Workflow for Ridgway’s Hawk Integrated Population Model
Brian W. Rolek
2024-11-20
Supplemental materials for: Rolek, B.W., McClure, CJW, Dunn, L., Curti, M., … Ridgway’s Hawk IPM and PVA
Contact information: rolek.brian@peregrinefund.org
Metadata, data, and scripts used in analyses can be found at https://github.com/The-Peregrine-Fund/XXXXX.
The full workflow below is visible as a html website at: https://the-peregrine-fund.github.io/XXXXX/.
A permanent archive and DOI is available at: https://zenodo.org/doi/XXXXX
1. Code for IPM and PVA
1.1 Submodels without Integration
library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
# Observations (po) = y
# 1 seen first-year (age=0, just before 1st b-day)
# 2 seen nonbreeder
# 3 seen breeder
# 4 not seen
# States (ps)
# 1 alive first-year
# 2 alive nonbreeder
# 3 alive breeder
# 4 dead
# Groups
# 1 wild-born
# 2 translocated and hacked
###################################################
# PARAMETERS
# phiFY: survival probability first year
# phiA: survival probability nonbreeders
# phiB: survival probability breeders
# psiFYB: recruitment probability from first-year to breeder
# psiAB: recruitment probability from nonbreeders to breeder
# psiBA: recruitment probability from breeder to nonbreeders
# pFY: resight probability first-years
# pA: resight probability nonbreeders
# pB: resight probability breeders
# lmu.prod: mean productivity (males and females) per territory on the log scale
# sds: standard deviations for multivariate normal random effects for time and site
# R: correlation coefficients for multivariate normal random effects for time and site
# lambda: population growth rate (derived)
# extinct: binary indicator of extirpation at a site
# gamma: coefficient of effect from nest treatments
# betas: coefficient of effect from translocations
# deltas: coefficient of effect from survey effort
# mus: overall means for survival, recruitment, and detections
# r and rr: "r" parameter for negative binomial distribution
# also described as omega in manuscript
#**********************
#* Model code
#**********************
mycode <- nimbleCode(
{
###################################################
# Priors and constraints
###################################################
# Survival, recruitment, and detection can be correlated
for (k in 1:8){
betas[k] ~ dunif(-20, 20) # prior for coefficients
} # k
for (j in 1:8){
for (s in 1:nsite){
lmus[j,s] <- logit(mus[j,s])
mus[j,s] ~ dbeta(1,1) # prior for means
}} # m population #s sex #h hacked
# # Temporal random effects and correlations among all sites, synchrony
# for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation estimated using the multivariate normal distribution
# R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
# Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
# U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
# # Multivariate normal for temporal variance
# for (t in 1:nyr){
# eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
# cholesky = U[1:p, 1:p], prec_param = 0)
# }
# Temporal random effects and correlations between sites
for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation estimated using the multivariate normal distribution
R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
# multivariate normal for temporal and site variance
for (t in 1:nyr){
for (s in 1:nsite){
eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
cholesky = U2[1:p2, 1:p2], prec_param = 0)
} } # s t
###############################
# Likelihood for productivity
###############################
# Priors for productivity
for (s in 1:nsite){
lmu.f[s] ~ dnorm(0, sd=5)
} # s
gamma ~ dunif(-20, 20)
rr ~ dexp(0.05)
# Productivity likelihood
for (k in 1:npairsobs){
f[k] ~ dnegbin(ppp[k], rr)
ppp[k] <- rr/(rr+mu.f[k])
log(mu.f[k]) <- lmu.f[site.pair[k]] +
gamma*treat.pair[k] +
#eps[9, year.pair[k] ] +
eta[9, site.pair[k], year.pair[k] ]
} # k
# Derive yearly brood size for population model
# Need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.pair is a matrix of indices for each site
for (t in 1:nyr){
for (s in 1:nsite){
for (xxx in 1:pair.end[t,s]){
fecmat[t,s,xxx] <- mu.f[ yrind.pair[xxx,t,s] ]
} # xxx
mn.f[t,s] <- mean( fecmat[t,s,1:pair.end[t,s]] )
}} # s t
# GOF for number of fledglings
for (k in 1:npairsobs){
f.obs[k] <- f[k] # observed counts
f.exp[k] <- mu.f[k] # expected counts adult breeder
f.rep[k] ~ dnegbin(ppp[k], rr) # expected counts
f.dssm.obs[k] <- abs( ( f.obs[k] - f.exp[k] ) / (f.obs[k]+0.001) )
f.dssm.rep[k] <- abs( ( f.rep[k] - f.exp[k] ) / (f.rep[k]+0.001) )
} # k
f.dmape.obs <- sum(f.dssm.obs[1:npairsobs])
f.dmape.rep <- sum(f.dssm.rep[1:npairsobs])
f.tvm.obs <- sd(brood[1:npairsobs])^2/mean(brood[1:npairsobs])
f.tvm.rep <- sd(f.rep[1:npairsobs])^2/mean(f.rep[1:npairsobs])
################################
# Likelihood for counts
################################
# Omitted. See IPM or PVA.
################################
# Likelihood for survival
################################
# Calculate yearly averages for sites for integration
for (t in 1:nyr){
for (s in 1:nsite){
for (xxxx in 1:surv.end[t,s]){
# Need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.surv is a matrix of indices for each site
phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
} # xxxx
mn.phiFY[ t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] )
mn.phiA[ t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
mn.phiB[ t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
mn.psiFYB[ t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
mn.psiAB[ t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
mn.psiBA[ t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
mn.pA[ t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
mn.pB[ t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
}}
for (i in 1:nind){
for (t in 1:nyr){
#Survival
logit(phiFY[i,t]) <- eta[1, site[i,t],t] + # eps[1,t] +
lmus[1, site[i,t]] + betas[1]*hacked[i] # first year
logit(phiA[i,t]) <- eta[2, site[i,t],t] +# eps[2,t] +
lmus[2, site[i,t]] + betas[2]*hacked[i] # nonbreeder
logit(phiB[i,t]) <- eta[3, site[i,t],t] +# eps[3,t] +
lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
#Recruitment
logit(psiFYB[i,t]) <- eta[4, site[i,t],t] +# eps[4,t] +
lmus[4, site[i,t]] + betas[4]*hacked[i] # first year to breeder
logit(psiAB[i,t]) <- eta[5, site[i,t],t] +# eps[5,t] +
lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
logit(psiBA[i,t]) <- #eta[6, site[i,t],t] + eps[6,t] +
lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
#Re-encounter
logit(pA[i,t]) <- eta[7, site[i,t],t] +# eps[7,t] +
lmus[7, site[i,t]] + betas[7]*hacked[i] # resight of nonbreeders
logit(pB[i,t]) <- eta[8, site[i,t],t] +# eps[8,t] +
lmus[8, site[i,t]] + betas[8]*hacked[i] # resight of breeders
}#t
}#i
# Define state-transition and observation matrices
for (i in 1:nind){
# Define probabilities of state S(t+1) given S(t)
for (t in first[i]:(nyr-1)){
ps[1,i,t,1] <- 0
ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t])
ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
ps[1,i,t,4] <- (1-phiFY[i,t])
ps[2,i,t,1] <- 0
ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
ps[2,i,t,4] <- (1-phiA[i,t])
ps[3,i,t,1] <- 0
ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
ps[3,i,t,4] <- (1-phiB[i,t])
ps[4,i,t,1] <- 0
ps[4,i,t,2] <- 0
ps[4,i,t,3] <- 0
ps[4,i,t,4] <- 1
# Define probabilities of O(t) given S(t)
po[1,i,t,1] <- 1
po[1,i,t,2] <- 0
po[1,i,t,3] <- 0
po[1,i,t,4] <- 0
po[2,i,t,1] <- 0
po[2,i,t,2] <- pA[i,t]
po[2,i,t,3] <- 0
po[2,i,t,4] <- (1-pA[i,t])
po[3,i,t,1] <- 0
po[3,i,t,2] <- 0
po[3,i,t,3] <- pB[i,t]
po[3,i,t,4] <- (1-pB[i,t])
po[4,i,t,1] <- 0
po[4,i,t,2] <- 0
po[4,i,t,3] <- 0
po[4,i,t,4] <- 1
} #t
} #i
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,first[i]] <- y[i,first[i]]
for (t in (first[i]+1):nyr){
# State process: draw S(t) given S(t-1)
z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
# Observation process: draw O(t) given S(t)
y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
} #t
} #i
} )
#**********************
#* Function to run model in NIMBLE
#**********************
run_ipm <- function(info, datl, constl, code){
library('nimble')
library('coda')
# helper function for multivariate normal
uppertri_mult_diag <- nimbleFunction(
run = function(mat = double(2), vec = double(1)) {
returnType(double(2))
p <- length(vec)
out <- matrix(nrow = p, ncol = p, init = FALSE)
for(i in 1:p){
out[ ,i] <- mat[ ,i] * vec[i]
}
return(out)
})
assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
params <- c(# pop growth
#"lambda",
# fecundity
"lmu.f", "gamma", "rr", "mn.f",
# survival
"mus", "lmus", "betas",
# abundance
# "NB", "NF", "NFY", "N",
# "r",
# "N", "Ntot",
# error terms
#"eps", "sds", "Ustar", "U", "R",
"eta", "sds2", "Ustar2", "U2", "R2",
# yearly summaries
'mn.phiFY', 'mn.phiA', 'mn.phiB',
'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
'mn.pA', 'mn.pB',
# goodness of fit
"f.dmape.obs", "f.dmape.rep",
"f.tvm.obs", "f.tvm.rep"
# "dmape.obs", "dmape.rep",
# "tvm.obs", "tvm.rep",
# "tturn.obs", "tturn.rep"
)
#n.chains=1; n.thin=200; n.iter=500000; n.burnin=300000
#n.chains=1; n.thin=50; n.iter=100000; n.burnin=50000
#n.chains=1; n.thin=10; n.iter=20000; n.burnin=10000
n.chains=1; n.thin=200; n.iter=500000; n.burnin=300000
mod <- nimbleModel(code,
constants = constl,
data = datl,
inits = info$inits,
buildDerivs = FALSE, # doesn't work when TRUE, no hope for HMC
calculate=T
)
cmod <- compileNimble(mod, showCompilerOutput = TRUE)
confhmc <- configureMCMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project = mod,
resetFunctions = TRUE,
showCompilerOutput = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T,
setSeed = info$seed)
return(post)
} # run_ipm function end
#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(10)
post <- parLapply(cl = this_cluster,
X = par_info,
fun = run_ipm,
datl = datl,
constl = constl,
code = mycode)
stopCluster(this_cluster)
save(post, mycode,
file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm_simp.rdata")
# save(post, mycode,
# file="C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_statespace.rdata")1.2 Integrated Population Model
library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
# Observations (po) = y
# 1 seen first-year (age=0, just before 1st b-day)
# 2 seen nonbreeder
# 3 seen breeder
# 4 not seen
# States (ps)
# 1 alive first-year
# 2 alive nonbreeder
# 3 alive breeder
# 4 dead
# Groups
# 1 wild-born
# 2 translocated and hacked
###################################################
# PARAMETERS
# phiFY: survival probability first year
# phiA: survival probability nonbreeders
# phiB: survival probability breeders
# psiFYB: recruitment probability from first-year to breeder
# psiAB: recruitment probability from nonbreeders to breeder
# psiBA: recruitment probability from breeder to nonbreeders
# pFY: resight probability first-years
# pA: resight probability nonbreeders
# pB: resight probability breeders
# lmu.prod: mean productivity (males and females) per territory on the log scale
# sds: standard deviations for multivariate normal random effects for time and site
# R: correlation coefficients for multivariate normal random effects for time and site
# lambda: population growth rate (derived)
# extinct: binary indicator of extirpation at a site
# gamma: coefficient of effect from nest treatments
# betas: coefficient of effect from translocations
# deltas: coefficient of effect from survey effort
# mus: overall means for survival, recruitment, and detections
# r and rr: "r" parameter for negative binomial distribution
# also described as omega in manuscript
#**********************
#* Model code
#**********************
mycode <- nimbleCode(
{
###################################################
# Priors and constraints
###################################################
# survival, recruitment, and detection can be correlated
for (k in 1:8){
betas[k] ~ dunif(-20, 20) # prior for translocations coefficients
deltas[k] ~ dunif(-20, 20) # prior for survey effort coefficients
} # k
for (j in 1:8){
for (s in 1:nsite){
lmus[j,s] <- logit(mus[j,s])
mus[j,s] ~ dbeta(1,1) # prior for overall means
}} #
# # Temporal random effects and correlations among all sites
# for (j in 1:p){ sds[j] ~ dexp(1) }# prior for temporal variation
# # estimated using the multivariate normal distribution
# R[1:p,1:p] <- t(Ustar[1:p,1:p]) %*% Ustar[1:p,1:p] # calculate rhos, correlation coefficients
# Ustar[1:p,1:p] ~ dlkj_corr_cholesky(eta=1.1, p=p) # Ustar is the Cholesky decomposition of the correlation matrix
# U[1:p,1:p] <- uppertri_mult_diag(Ustar[1:p, 1:p], sds[1:p])
# # multivariate normal for temporal variance
# for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
# eps[1:p,t] ~ dmnorm(mu.zeroes[1:p],
# cholesky = U[1:p, 1:p], prec_param = 0)
# }
# Temporal random effects and correlations between sites
for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal variation
# estimated using the multivariate normal distribution
R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2])
# multivariate normal for temporal variance
for (t in 1:nyr){ # survival params only have nyr-1, no problem to simulate from however
for (s in 1:nsite){
eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
cholesky = U2[1:p2, 1:p2], prec_param = 0)
} } # s t
#######################
# Derived params
#######################
for (s in 1:nsite){
for (t in 1:(nyr-1)){
lambda[t, s] <- Ntot[t+1, s]/(Ntot[t, s])
loglambda[t, s] <- log(lambda[t, s])
}} #t
###############################
# Likelihood for productivity
###############################
# Priors
for (s in 1:nsite){
lmu.prod[s] ~ dnorm(0, sd=5)
} # s
gamma ~ dunif(-20, 20)
rr ~ dexp(0.05)
# Productivity likelihood
for (k in 1:npairsobs){
prod[k] ~ dnegbin(ppp[k], rr)
ppp[k] <- rr/(rr+mu.prod[k])
log(mu.prod[k]) <- lmu.prod[site.pair[k]] +
gamma*treat.pair[k] +
#eps[9, year.pair[k] ] +
eta[9, site.pair[k], year.pair[k] ]
} # k
# Derive yearly productivity for population model
# need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.pair is a matrix of indices for each site
for (t in 1:nyr){
for (s in 1:nsite){
for (xxx in 1:pair.end[t,s]){
prodmat[t,s,xxx] <- mu.prod[ yrind.pair[xxx,t,s] ]
} # xxx
mn.prod[t,s] <- mean( prodmat[t,s,1:pair.end[t,s]] )
}} # s t
# GOF for productivity
for (k in 1:npairsobs){
f.obs[k] <- prod[k] # observed counts
f.exp[k] <- mu.prod[k] # expected counts adult breeder
f.rep[k] ~ dnegbin(ppp[k], rr) # expected counts
f.dssm.obs[k] <- abs( ( f.obs[k] - f.exp[k] ) / (f.obs[k]+0.001) )
f.dssm.rep[k] <- abs( ( f.rep[k] - f.exp[k] ) / (f.rep[k]+0.001) )
} # k
f.dmape.obs <- sum(f.dssm.obs[1:npairsobs])
f.dmape.rep <- sum(f.dssm.rep[1:npairsobs])
f.tvm.obs <- sd(brood[1:npairsobs])^2/mean(brood[1:npairsobs])
f.tvm.rep <- sd(f.rep[1:npairsobs])^2/mean(f.rep[1:npairsobs])
################################
# Likelihood for counts
################################
# Abundance for year=1
for (v in 1:7){
for (s in 1:nsite){
# subtract one to allow dcat to include zero
N[v, 1, s] <- N2[v, 1, s] - 1
N2[v, 1, s] ~ dcat(pPrior[v, 1:s.end[v,s], s]) # Priors differ for FYs and adults
}} # s t
# Abundance for years > 1
for (t in 1:(nyr-1)){
for (s in 1:nsite){
# Number of wild born juvs
N[1, t+1, s] ~ dpois( (NFY[t, s]*mn.phiFY[t, s]*mn.psiFYB[t, s] + # first year breeders
NF[t, s]*mn.phiA[t, s]*mn.psiAB[t, s] + # nonbreeders to breeders
NB[t, s]*mn.phiB[t, s]*(1-mn.psiBA[t, s])) # breeders remaining
*mn.prod[t+1, s]/2 ) # end Poisson
# Abundance of nonbreeders
N[2, t+1, s] ~ dbin(mn.phiFY[t, s]*(1-mn.psiFYB[t, s]), NFY[t, s]) # Nestlings to nonbreeders
N[3, t+1, s] ~ dbin(mn.phiA[t, s]*(1-mn.psiAB[t, s]), NF[t, s]) # Nonbreeders to nonbreeders
N[4, t+1, s] ~ dbin(mn.phiB[t, s]*mn.psiBA[t, s], NB[t, s]) # Breeders to nonbreeders
# Abundance of breeders
N[5, t+1, s] ~ dbin(mn.phiFY[t, s]*mn.psiFYB[t, s], NFY[t, s]) # Nestlings to second year breeders
N[6, t+1, s] ~ dbin(mn.phiA[t, s]*mn.psiAB[t, s], NF[t, s]) # Nonbreeder to breeder
N[7, t+1, s] ~ dbin(mn.phiB[t, s]*(1-mn.psiBA[t, s]), NB[t, s]) # Breeder to breeder
}} # s t
# Count likelihoods, state-space model, and observation process
for (t in 1:nyr){
for (s in 1:nsite){
countsAdults[t, s] ~ dpois(lamAD[t, s]) # adult females
constraint_data[t, s] ~ dconstraint( (N[1, t, s] + hacked.counts[t, s]) >= 0 ) # constrain N1+hacked.counts to be >=0
NFY[t, s] <- N[1, t, s] + hacked.counts[t, s] # Transfers translocated first-year females
NF[t, s] <- sum(N[2:4, t, s]) # number of adult nonbreeder females
NB[t, s] <- sum(N[5:7, t, s]) # number of adult breeder females
NAD[t, s] <- NF[t, s] + NB[t, s] # number of adults
Ntot[t, s] <- sum(N[1:7, t, s]) # total number of females
log(lamAD[t, s]) <- log(NAD[t, s]) + deltas[1]*effort2[t, s] + deltas[2]*effort2[t, s]^2
log(lamFY[t, s]) <- log(N[1, t, s]) + deltas[3]*effort2[t, s] + deltas[4]*effort2[t, s]^2
}# s
# First-years at different sites have different distributions
# for better model fit
countsFY[t, 1] ~ dpois(lamFY[t, 1]) # doesn't have any zeroes so poisson
countsFY[t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
pp[t] <- r/(r+(lamFY[t, 2] ))
} # t
r ~ dexp(0.05)
###################
# Assess GOF of the state-space models for counts
# Step 1: Compute statistic for observed data
# Step 2: Use discrepancy measure: mean absolute error
# Step 3: Use test statistic: number of turns
###################
for (t in 1:nyr){
c.repFY[t, 1] ~ dpois( lamFY[t, 1] )
c.repFY[t, 2] ~ dnegbin( pp[t], r )
for (s in 1:nsite){
c.expAD[t, s] <- lamAD[t, s] # expected counts adult breeder
c.expFY[t, s] <- lamFY[t, s]
c.obsAD[t, s] <- countsAdults[t, s]
c.obsFY[t, s] <- countsFY[t, s] # first year
c.repAD[t, s] ~ dpois( lamAD[t, s] ) # simulated counts
dssm.obsAD[t, s] <- abs( ( (c.obsAD[t, s]) - (c.expAD[t, s]) ) / (c.obsAD[t, s]+0.001) )
dssm.obsFY[t, s] <- abs( ( (c.obsFY[t, s]) - (c.expFY[t, s]) ) / (c.obsFY[t, s]+0.001) )
dssm.repAD[t, s] <- abs( ( (c.repAD[t, s]) - (c.expAD[t, s]) ) / (c.repAD[t, s]+0.001) )
dssm.repFY[t, s] <- abs( ( (c.repFY[t, s]) - (c.expFY[t, s]) ) / (c.repFY[t, s]+0.001) )
}} # t
dmape.obs[1] <- sum(dssm.obsAD[1:nyr, 1:nsite])
dmape.obs[2] <- sum(dssm.obsFY[1:nyr, 1:nsite])
dmape.rep[1] <- sum(dssm.repAD[1:nyr, 1:nsite])
dmape.rep[2] <- sum(dssm.repFY[1:nyr, 1:nsite])
tvm.obs[1] <- sd(c.obsB[1:nyr, 1:nsite])^2/mean(c.obsB[1:nyr, 1:nsite])
tvm.obs[2] <- sd(c.obsFY[1:nyr, 1:nsite])^2/mean(c.obsFY[1:nyr, 1:nsite])
tvm.rep[1] <- sd(c.repB[1:nyr, 1:nsite])^2/mean(c.repB[1:nyr, 1:nsite])
tvm.rep[2] <- sd(c.repFY[1:nyr, 1:nsite])^2/mean(c.repFY[1:nyr, 1:nsite])
# # Test statistic for number of turns
for (s in 1:nsite){
for (t in 1:(nyr-2)){
tt1.obsAD[t, s] <- step(c.obsAD[t+2, s] - c.obsAD[t+1, s])
tt2.obsAD[t, s] <- step(c.obsAD[t+1, s] - c.obsAD[t, s])
tt3.obsAD[t, s] <- equals(tt1.obsAD[t, s] + tt2.obsAD[t, s], 1)
tt1.obsFY[t, s] <- step(c.obsFY[t+2, s] - c.obsFY[t+1, s])
tt2.obsFY[t, s] <- step(c.obsFY[t+1, s] - c.obsFY[t, s])
tt3.obsFY[t, s] <- equals(tt1.obsFY[t, s] + tt2.obsFY[t, s], 1)
}} # t
tturn.obs[1] <- sum(tt3.obsAD[1:(nyr-2), 1:nsite])
tturn.obs[2] <- sum(tt3.obsFY[1:(nyr-2), 1:nsite])
for (s in 1:nsite){
for (t in 1:(nyr-2)){
tt1.repAD[t, s] <- step(c.repAD[t+2, s] - c.repAD[t+1, s])
tt2.repAD[t, s] <- step(c.repAD[t+1, s] - c.repAD[t, s])
tt3.repAD[t, s] <- equals(tt1.repAD[t, s] + tt2.repAD[t, s], 1)
tt1.repFY[t, s] <- step(c.repFY[t+2, s] - c.repFY[t+1, s])
tt2.repFY[t, s] <- step(c.repFY[t+1, s] - c.repFY[t, s])
tt3.repFY[t, s] <- equals(tt1.repFY[t, s] + tt2.repFY[t, s], 1)
}} # t
tturn.rep[1] <- sum(tt3.repAD[1:(nyr-2), 1:nsite])
tturn.rep[2] <- sum(tt3.repFY[1:(nyr-2), 1:nsite])
################################
# Likelihood for survival
################################
# Calculate averages for sites each year for integration
for (t in 1:nyr){
for (s in 1:nsite){
for (xxxx in 1:surv.end[t,s]){
# Reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.surv is a matrix of indices for each site
phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
} # xxxx
mn.phiFY[ t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] )
mn.phiA[ t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
mn.phiB[ t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
mn.psiFYB[ t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
mn.psiAB[ t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
mn.psiBA[ t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
mn.pA[ t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
mn.pB[ t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
}}
for (i in 1:nind){
for (t in 1:nyr){
#Survival
logit(phiFY[i,t]) <- eta[1, site[i,t],t] + # eps[1,t] +
lmus[1, site[i,t]] + betas[1]*hacked[i] # first year
logit(phiA[i,t]) <- eta[2, site[i,t],t] + #eps[2,t] +
lmus[2, site[i,t]] + betas[2]*hacked[i] # nonbreeder
logit(phiB[i,t]) <- eta[3, site[i,t],t] + #eps[3,t] +
lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
#Recruitment
logit(psiFYB[i,t]) <- eta[4, site[i,t],t] + #eps[4,t] +
lmus[4, site[i,t]] + betas[4]*hacked[i] # first year to breeder
logit(psiAB[i,t]) <- eta[5, site[i,t],t] + #eps[5,t] +
lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
logit(psiBA[i,t]) <- #eta[6, site[i,t],t] + eps[6,t] +
lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
#Re-encounter
logit(pA[i,t]) <- eta[7, site[i,t],t] + #eps[7,t] +
lmus[7, site[i,t]] + betas[7]*hacked[i] +
deltas[5]*effort2[t, site[i,t]] + deltas[6]*effort2[t, site[i,t]]^2# resight of nonbreeders
logit(pB[i,t]) <- eta[8, site[i,t],t] + #eps[8,t] +
lmus[8, site[i,t]] + betas[8]*hacked[i] +
deltas[7]*effort2[t, site[i,t]] + deltas[8]*effort2[t, site[i,t]]^2# resight of breeders
}#t
}#i
# Define state-transition and observation matrices
for (i in 1:nind){
# Define probabilities of state S(t+1) given S(t)
for (t in first[i]:(nyr-1)){
ps[1,i,t,1] <- 0
ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t])
ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
ps[1,i,t,4] <- (1-phiFY[i,t])
ps[2,i,t,1] <- 0
ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
ps[2,i,t,4] <- (1-phiA[i,t])
ps[3,i,t,1] <- 0
ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
ps[3,i,t,4] <- (1-phiB[i,t])
ps[4,i,t,1] <- 0
ps[4,i,t,2] <- 0
ps[4,i,t,3] <- 0
ps[4,i,t,4] <- 1
# Define probabilities of O(t) given S(t)
po[1,i,t,1] <- 1
po[1,i,t,2] <- 0
po[1,i,t,3] <- 0
po[1,i,t,4] <- 0
po[2,i,t,1] <- 0
po[2,i,t,2] <- pA[i,t]
po[2,i,t,3] <- 0
po[2,i,t,4] <- (1-pA[i,t])
po[3,i,t,1] <- 0
po[3,i,t,2] <- 0
po[3,i,t,3] <- pB[i,t]
po[3,i,t,4] <- (1-pB[i,t])
po[4,i,t,1] <- 0
po[4,i,t,2] <- 0
po[4,i,t,3] <- 0
po[4,i,t,4] <- 1
} #t
} #i
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,first[i]] <- y[i,first[i]]
for (t in (first[i]+1):nyr){
# State process: draw S(t) given S(t-1)
z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
# Observation process: draw O(t) given S(t)
y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
} #t
} #i
} )
#**********************
#* Function to run model in NIMBLE
#**********************
run_ipm <- function(info, datl, constl, code){
library('nimble')
library('coda')
# helper function for multivariate normal
uppertri_mult_diag <- nimbleFunction(
run = function(mat = double(2), vec = double(1)) {
returnType(double(2))
p <- length(vec)
out <- matrix(nrow = p, ncol = p, init = FALSE)
for(i in 1:p){
out[ ,i] <- mat[ ,i] * vec[i]
}
return(out)
})
assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
params <- c(
# pop growth
"lambda",
# fecundity
"lmu.prod", "gamma", "rr", "mn.prod",
# survival
"mus", "lmus", "betas", "deltas",
# abundance
"NB", "NF", "NFY", "N", "NAD",
"r",
"N", "Ntot",
# error terms
#"eps", "sds", "Ustar", "U", "R",
"eta", "sds2", "Ustar2", "U2", "R2",
# yearly summaries
'mn.phiFY', 'mn.phiA', 'mn.phiB',
'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
'mn.pA', 'mn.pB',
# goodness of fit
"f.dmape.obs", "f.dmape.rep",
"f.tvm.obs", "f.tvm.rep",
"dmape.obs", "dmape.rep",
"tvm.obs", "tvm.rep",
"tturn.obs", "tturn.rep"
)
#n.chains=1; n.thin=1; n.iter=500; n.burnin=100
#n.chains=1; n.thin=200; n.iter=800000; n.burnin=600000
n.chains=1; n.thin=200; n.iter=600000; n.burnin=400000
mod <- nimbleModel(code,
constants = constl,
data = datl,
inits = info$inits,
buildDerivs = FALSE, # doesn't work when TRUE, no hope for HMC
calculate=T
)
cmod <- compileNimble(mod, showCompilerOutput = TRUE)
confhmc <- configureMCMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project = mod,
resetFunctions = TRUE,
showCompilerOutput = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T,
setSeed = info$seed)
return(post)
} # run_ipm function end
#*****************
#* Run chains in parallel
#*****************
par_info <- par_info[1:10]
this_cluster <- makeCluster(10)
post <- parLapply(cl = this_cluster,
X = par_info,
fun = run_ipm,
datl = datl,
constl = constl,
code = mycode)
stopCluster(this_cluster)
# save(post, mycode,
# file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm_shortrun.rdata")
save(post, mycode,
file="/bsuscratch/brianrolek/riha_ipm/outputs/ipm_longrun.rdata")1.3 Population Viability Analysis
library('nimble')
library('parallel')
library ('coda')
load("/bsuscratch/brianrolek/riha_ipm/data.rdata")
source("/bsuscratch/brianrolek/riha_ipm/MCMCvis.R")
load("/bsuscratch/brianrolek/riha_ipm/outputs/ipm_longrun.rdata")
# load("C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//outputs//ipm_longrun.rdata")
# load("data//data.rdata")
# library ("MCMCvis")
#**********************
#* Set data to reflect PVA scenarios
#**********************
constl$K <- 50 # number of future years
cpus <- 10 # number of processors
constl$SC <- 36 # number of PVA scenarios
datl$constraint_data <- rbind(datl$constraint_data, array(1, dim=c(constl$K,2)) ) # to help constrain FYs born + hacked to be positive
constl$effort2 <- rbind(constl$effort2, array(0, dim=c(constl$K,2), dimnames=list(2024:(2024+constl$K-1), c("LH", "PC"))))
constl$num.treated <- rep( c(0, 15, 30, 45, 0, 15, 30, 45, 0, 15, 30, 45), each=3 ) # 100s sub in for All but are over-ridden in model code
num.hacked <- rep( c(0, 0, 0, 0, 5, 5, 5, 5, 10, 10, 10, 10), each=3 )
surv.diff <- array(NA, dim=c(constl$SC, 3, 2), dimnames=list(scenario=1:constl$SC, stage=c('FY', 'NB', 'B'), site=c('LH', 'PC')))
surv.diff[,1,] <- matrix(c(rep( c(0, 0.130, 0.260), 12),
rep( c(0, 0.146, 0.292), 12)), ncol=2)
surv.diff[,2,] <- matrix(c(rep( c(0, 0.006, 0.012), 12),
rep( c(0, 0.022, 0.044), 12)), ncol=2)
surv.diff[,3,] <- matrix(c(rep( c(0, 0.016, 0.032), 12),
rep( c(0, 0.026, 0.052), 12)), ncol=2)
constl$surv.diff <- surv.diff
hacked.counts <- array(0, dim=c(constl$SC, constl$nyr+constl$K, 2))
for (sc in 1:constl$SC){
hacked.counts[sc,1:constl$nyr, 1:2] <- constl$hacked.counts[, 1:2]
hacked.counts[sc,14:(constl$nyr+constl$K), 1] <- -num.hacked[sc]
hacked.counts[sc,14:(constl$nyr+constl$K), 2] <- 0
}
constl$hacked.counts <- hacked.counts
#**********************
#* Specify priors taken from posterior of IPM
#**********************
#* This helps ensure that all chains run
# Identify chains with NAs that failed to initialize
out <- lapply(post, as.mcmc)
NAlist <- c()
for (i in 1:length(out)){
NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
out <- out[!NAlist]
outp <- MCMCpstr(out, type="chains")
!NAlist
Ni.func <- function (){
# randomly select inits from ones that work
# Ni <- outp$N[1:7,1:constl$nyr,1:2,
# sample(seq(1, 10000, by=100), 1, replace = F)]
# take the means of the posterior for inits
Ni <- apply(outp$N, c(1,2,3), mean) |> round()
Ni.pva <- array(NA, dim=c(constl$SC, dim(Ni)+c(0,constl$K,0)))
for (sc in 1:constl$SC){
Ni.pva[sc, 1:7, 1:constl$nyr, 1:2] <- Ni
for (t in constl$nyr:(constl$nyr+constl$K)){
for (s in 1:2){
Ni.pva[sc, 1:7, t, s] <- Ni[1:7, 13, s]
}}} # t sc
return(Ni.pva)
} # function
u2 <- apply(outp$Ustar2, c(1,2), mean)
inits.func.pva <- function (){
list(
# fecundity
lmu.prod = apply(outp$lmu.prod, 1, mean),
gamma = mean(outp$gamma),
rr = mean(outp$rr),
# survival
z = par_info[[1]]$inits$z,
mus = apply(outp$mus, c(1,2), mean),
betas = apply(outp$betas, 1, mean),
deltas = apply(outp$deltas, 1, mean),
sds2 = apply(outp$sds2, 1, mean),
Ustar2 = u2,
# counts
countsAdults= matrix(c(374, 335, 305, 295, rep(NA, length(2015:2023)), rep(NA, length(2011:2023)) ), nrow=13),
r = mean(outp$r),
N = Ni.func()
)}
# Set seed then save for reproducibility
# then randomly draw for each chain
set.seed(1)
seeds <- sample(1:1000000, size=cpus, replace=FALSE)
par_info_pva <- list()
for (i in 1:cpus){
par_info_pva[[i]] <- list(seed=seeds[i], inits = inits.func.pva())
}
#**********************
#* Parameter descriptions
#**********************
# Mark-resight-recovery data
# Observations (po) = y
# 1 seen first-year (age=0, just before 1st b-day)
# 2 seen nonbreeder
# 3 seen breeder
# 4 not seen
# States (ps)
# 1 alive first-year
# 2 alive nonbreeder
# 3 alive breeder
# 4 dead
# Groups
# 1 wild-born
# 2 translocated and hacked
###################################################
# PARAMETERS
# phiFY: survival probability first year
# phiA: survival probability nonbreeders
# phiB: survival probability breeders
# psiFYB: recruitment probability from first-year to breeder
# psiAB: recruitment probability from nonbreeders to breeder
# psiBA: recruitment probability from breeder to nonbreeders
# pFY: resight probability first-years
# pA: resight probability nonbreeders
# pB: resight probability breeders
# lmu.prod: mean productivity (males and females) per territory on the log scale
# sds: standard deviations for multivariate normal random effects for time and site
# R: correlation coefficients for multivariate normal random effects for time and site
# lambda: population growth rate (derived)
# extinct: binary indicator of extirpation at a site
# gamma: coefficient of effect from nest treatments
# betas: coefficient of effect from translocations
# deltas: coefficient of effect from survey effort
# mus: overall means for survival, recruitment, and detections
# r and rr: "r" parameter for negative binomial distribution
# also described as omega in manuscript
#**********************
#* Model code
#**********************
mycode <- nimbleCode(
{
###################################################
# Priors and constraints
###################################################
# survival, recruitment, and detection can be correlated
for (k in 1:8){
betas[k] ~ dunif(-20, 20) # prior for coefficients
deltas[k] ~ dunif(-20, 20)
} # k
for (j in 1:8){
for (s in 1:nsite){
lmus[j,s] <- logit(mus[j,s])
mus[j,s] ~ dbeta(1,1) # prior for means
}} # m population #s sex #h hacked
# Temporal and site random effects using the multivariate normal distribution
for (jj in 1:p2){ sds2[jj] ~ dexp(1) }# prior for temporal and site variation
R2[1:p2,1:p2] <- t(Ustar2[1:p2,1:p2]) %*% Ustar2[1:p2,1:p2] # calculate rhos, correlation coefficients
Ustar2[1:p2,1:p2] ~ dlkj_corr_cholesky(eta=1.1, p=p2) # Ustar is the Cholesky decomposition of the correlation matrix
U2[1:p2,1:p2] <- uppertri_mult_diag(Ustar2[1:p2, 1:p2], sds2[1:p2]) # custom function taken from nimble listserv (specified below)
# multivariate normal for temporal variance
for (t in 1:(nyr+K) ){
for (s in 1:nsite){
eta[1:p2,s,t] ~ dmnorm(mu.zeroes2[1:p2],
cholesky = U2[1:p2, 1:p2], prec_param = 0)
} } # s t
#######################
# Derived params
#######################
for(sc in 1:SC){
for (s in 1:nsite){
for (t in 1:(nyr-1+K)){
lambda[sc, t, s] <- Ntot[sc, t+1, s]/(Ntot[sc, t, s]) # population groewth rate
loglambda[sc, t, s] <- log(lambda[sc, t, s]) # r
}} #t
for (s in 1:nsite){
for (t in 1:K){
# quasi-extinction probabilities given population segments, N total, adults, and breeders
extinct[sc, t, s] <- equals(Ntot[sc, nyr+t, s], 0)
extinctAD[sc, t, s] <- equals(NAD[sc, nyr+t, s], 0)
extinctB[sc, t, s] <- equals(NB[sc, nyr+t, s], 0)
}}
} # sc
###############################
# Likelihood for productivity
###############################
# Priors
for (s in 1:nsite){
lmu.prod[s] ~ dnorm(0, sd=5)
} # s
gamma ~ dunif(-20, 20)
rr ~ dexp(0.05)
# Productivity likelihood
for (k in 1:npairsobs){
prod[k] ~ dnegbin(ppp[k], rr)
ppp[k] <- rr/(rr+mu.prod[k])
log(mu.prod[k]) <- lmu.prod[site.pair[k]] +
gamma*treat.pair[k] +
eta[9, site.pair[k], year.pair[k] ]
} # k
# Derive yearly productivity for population model
# need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.pair is a matrix of indices for each site
for (t in 1:nyr){
for (s in 1:nsite){
for (xxx in 1:pair.end[t,s]){
prodmat[t,s,xxx] <- mu.prod[ yrind.pair[xxx,t,s] ]
} # xxx
mn.prod[1,t,s] <- mean( prodmat[t,s,1:pair.end[t,s]] )
for(sc in 2:SC){
mn.prod[sc, t, s] <- mn.prod[1,t,s]
}}} # s t
# Future projections- fecundity
for(sc in 1:SC){
for (t in (nyr+1):(nyr+K) ){
for (s in 1:nsite){
perc.treat[sc, t, s] <- #
step( NB[sc, t, s]-num.treated[sc] ) * (num.treated[sc] + 0.001)/(NB[sc, t, s]+0.001) + # num.treated > NB, calculate proportion
1-step( NB[sc, t, s]-num.treated[sc] ) # num.treated < NB set to 1
log(mn.prod[sc, t, s]) <- lmu.prod[s] +
gamma*perc.treat[sc, t, s] + # treat.pair set to one here.
eta[9, s, t]
}} } # s t sc
################################
# Likelihood for counts
################################
# Abundance for year=1
for (v in 1:7){
for (s in 1:nsite){
# Subtract one to allow dcat to include zero
N[1, v, 1, s] <- N2[1, v, 1, s] - 1
N2[1, v, 1, s] ~ dcat(pPrior[v, 1:s.end[v,s], s])
# Assign estimates for years with data to all scenarios
for (sc in 2:SC){
N[sc, v, 1, s] <- N[1, v, 1, s]
}
}} # s t
# Abundance for years > 1
for (t in 1:(nyr-1)){
for (s in 1:nsite){
# Number of wild born juvs
N[1, 1, t+1, s] ~ dpois( (NFY[1, t, s]*mn.phiFY[1, t, s]*mn.psiFYB[1, t, s] + # first year breeders
NF[1, t, s]*mn.phiA[1, t, s]*mn.psiAB[1, t, s] + # nonbreeders to breeders
NB[1, t, s]*mn.phiB[1, t, s]*(1-mn.psiBA[1, t, s])) # breeders remaining
*(mn.prod[1, t+1, s]/2) ) # end Poisson
# Abundance of nonbreeders
N[1, 2, t+1, s] ~ dbin(mn.phiFY[1, t, s]*(1-mn.psiFYB[1, t, s]), NFY[1, t, s]) # Nestlings to second year nonbreeders
N[1, 3, t+1, s] ~ dbin(mn.phiA[1, t, s]*(1-mn.psiAB[1, t, s]), NF[1, t, s]) # Nonbreeders to nonbreeders
N[1, 4, t+1, s] ~ dbin(mn.phiB[1, t, s]*mn.psiBA[1, t, s], NB[1, t, s]) # Breeders to nonbreeders
# Abundance of breeders
N[1, 5, t+1, s] ~ dbin(mn.phiFY[1, t, s]*mn.psiFYB[1, t, s], NFY[1, t, s]) # Nestlings to second year breeders
N[1, 6, t+1, s] ~ dbin(mn.phiA[1, t, s]*mn.psiAB[1, t, s], NF[1, t, s]) # Nonbreeder to breeder
N[1, 7, t+1, s] ~ dbin(mn.phiB[1, t, s]*(1-mn.psiBA[1, t, s]), NB[1, t, s]) # Breeder to breeder
# Assign estimates for years with data to
# all scenarios
for (v in 1:7){
for (sc in 2:SC){
N[sc, v, t+1, s] <- N[1, v, t+1, s]
} } # sc
}} # s t
# Future projections- population model
for (sc in 1:SC){
for (t in nyr:(nyr+K-1)){
for (s in 1:nsite){
# Number of wild born juvs
N[sc, 1, t+1, s] ~ dpois( (NFY[sc, t, s]*mn.phiFY[sc, t, s]*mn.psiFYB[sc, t, s] + # first year breeders
NF[sc, t, s]*mn.phiA[sc, t, s]*mn.psiAB[sc, t, s] + # nonbreeders to breeders
NB[sc, t, s]*mn.phiB[sc, t, s]*(1-mn.psiBA[sc, t, s])) # breeders remaining
*(mn.prod[sc, t+1, s]/2) ) # end Poisson
# Abundance of nonbreeders
## Second year nonbreeders
N[sc, 2, t+1, s] ~ dbin(mn.phiFY[sc, t, s]*(1-mn.psiFYB[sc, t, s]), NFY[sc, t, s]) # Nestlings to second year nonbreeders
## Adult nonbreeders
N[sc, 3, t+1, s] ~ dbin(mn.phiA[sc, t, s]*(1-mn.psiAB[sc, t, s]), NF[sc, t, s]) # Nonbreeders to nonbreeders
N[sc, 4, t+1, s] ~ dbin(mn.phiB[sc, t, s]*mn.psiBA[sc, t, s], NB[sc, t, s]) # Breeders to nonbreeders
# Abundance of breeders
## Second year breeders
N[sc, 5, t+1, s] ~ dbin(mn.phiFY[sc, t, s]*mn.psiFYB[sc, t, s], NFY[sc, t, s]) # Nestlings to second year breeders
## Adult breeders
N[sc, 6, t+1, s] ~ dbin(mn.phiA[sc, t, s]*mn.psiAB[sc, t, s], NF[sc, t, s]) # Nonbreeder to breeder
N[sc, 7, t+1, s] ~ dbin(mn.phiB[sc, t, s]*(1-mn.psiBA[sc, t, s]), NB[sc, t, s]) # Breeder to breeder
}} # s t
} # sc
# Count likelihoods, state-space model, and observation process
for (t in 1:nyr){
for (s in 1:nsite){
constraint_data[t, s] ~ dconstraint( (N[1, 1, t, s] + hacked.counts[1, t, s]) >= 0 ) # Transfers translocated first-year females
NFY[1, t, s] <- N[1, 1, t, s] + hacked.counts[1, t, s] # Transfers translocated first-year females
NF[1, t, s] <- sum(N[1, 2:4, t, s]) # number of adult nonbreeder females
NB[1, t, s] <- sum(N[1, 5:7, t, s]) # number of adult breeder females
NAD[1, t, s] <- NF[1, t, s] + NB[1, t, s] # number of adults
Ntot[1, t, s] <- sum(N[1, 1:7, t, s]) # total number of females
countsAdults[t, s] ~ dpois(lamAD[1, t, s]) # adult females
# constrain N1+hacked.counts to be >=0
log(lamAD[1, t, s]) <- log(NAD[1, t, s]) + deltas[1]*effort2[t, s] + deltas[2]*effort2[t, s]^2
log(lamFY[1, t, s]) <- log(N[1, 1, t, s]) + deltas[3]*effort2[t, s] + deltas[4]*effort2[t, s]^2
for (sc in 2:SC){
lamAD[sc, t, s] <- lamAD[1, t, s]
lamFY[sc, t, s] <- lamFY[1, t, s]
NFY[sc, t, s] <- NFY[1, t, s]
NF[sc, t, s] <- NF[1, t, s]
NB[sc, t, s] <- NB[1, t, s]
NAD[sc, t, s] <- NAD[1, t, s]
Ntot[sc, t, s] <- Ntot[1, t, s]
} # sc
}# s
# First-years at different sites have different distributions
# for better model fit
countsFY[t, 1] ~ dpois(lamFY[1, t, 1]) # doesn't have any zeroes so poisson
countsFY[t, 2] ~ dnegbin(pp[t], r) # first year females, includes translocated/hacked
pp[t] <- r/(r+(lamFY[1, t, 2] ))
} # t
r ~ dexp(0.05)
# Future projections of counts
for(sc in 1:SC){
for (t in (nyr+1):(nyr+K)){
for (s in 1:nsite){
# constrain N1+hacked.counts to be >=0, and allow it to shrink as the population shrinks
NFY[sc, t, s] <- step( (N[sc, 1, t, s] + hacked.counts[sc, t, s]) ) * (N[sc, 1, t, s] + hacked.counts[sc, t, s])
NF[sc, t, s] <- sum(N[sc, 2:4, t, s]) # number of adult nonbreeder females
NB[sc, t, s] <- sum(N[sc, 5:7, t, s]) # number of adult breeder females
NAD[sc, t, s] <- NF[sc, t, s] + NB[sc, t, s] # number of adults
Ntot[sc, t, s] <- sum(N[sc, 1:7, t, s]) # total number of females
}# s
} }# t sc
################################
# Likelihood for survival
################################
# Calculate yearly averages for integration
for (t in 1:(nyr-1)){
for (s in 1:nsite){
for (xxxx in 1:surv.end[t,s]){
# need to reorder because nimble doesn't
# handle nonconsecutive indices
# yrind.surv is a matrix of indices for each site
phiFY2[ t, s, xxxx] <- phiFY[ yrind.surv[xxxx,t,s], t]
phiA2[ t, s, xxxx] <- phiA[ yrind.surv[xxxx,t,s], t]
phiB2[ t, s, xxxx] <- phiB[ yrind.surv[xxxx,t,s], t]
psiFYB2[ t, s, xxxx] <- psiFYB[ yrind.surv[xxxx,t,s], t]
psiAB2[ t, s, xxxx] <- psiAB[ yrind.surv[xxxx,t,s], t]
psiBA2[ t, s, xxxx] <- psiBA[ yrind.surv[xxxx,t,s], t]
pA2[ t, s, xxxx] <- pA[ yrind.surv[xxxx,t,s], t]
pB2[ t, s, xxxx] <- pB[ yrind.surv[xxxx,t,s], t]
} # xxxx
mn.phiFY[1, t, s] <- mean( phiFY2[ t, s, 1:surv.end[t,s] ] )
mn.phiA[1, t, s] <- mean( phiA2[ t, s, 1:surv.end[t,s] ] )
mn.phiB[1, t, s] <- mean( phiB2[ t, s, 1:surv.end[t,s] ] )
mn.psiFYB[1, t, s] <- mean( psiFYB2[ t, s, 1:surv.end[t,s] ] )
mn.psiAB[1, t, s] <- mean( psiAB2[ t, s, 1:surv.end[t,s] ] )
mn.psiBA[1, t, s] <- mean( psiBA2[ t, s, 1:surv.end[t,s] ] )
mn.pA[1, t, s] <- mean( pA2[ t, s, 1:surv.end[t,s] ] )
mn.pB[1, t, s] <- mean( pB2[ t, s, 1:surv.end[t,s] ] )
for (sc in 2:SC){
mn.phiFY[sc, t, s] <- mn.phiFY[1, t, s]
mn.phiA[sc, t, s] <- mn.phiA[1, t, s]
mn.phiB[sc, t, s] <- mn.phiB[1, t, s]
mn.psiFYB[sc, t, s] <- mn.psiFYB[1, t, s]
mn.psiAB[sc, t, s] <- mn.psiAB[1, t, s]
mn.psiBA[sc, t, s] <- mn.psiBA[1, t, s]
mn.pA[sc, t, s] <- mn.pA[1, t, s]
mn.pB[sc, t, s] <- mn.pB[1, t, s]
}
}}
for (i in 1:nind){
for (t in 1:(nyr-1)){
#Survival
logit(phiFY[i,t]) <- eta[1, site[i,t],t] +
lmus[1, site[i,t]] + betas[1]*hacked[i] # first year
logit(phiA[i,t]) <- eta[2, site[i,t],t] +
lmus[2, site[i,t]] #+ betas[2]*hacked[i] # nonbreeder
logit(phiB[i,t]) <- eta[3, site[i,t],t] +
lmus[3, site[i,t]] + betas[3]*hacked[i] # breeder
#Recruitment
logit(psiFYB[i,t]) <- eta[4, site[i,t],t] +
lmus[4, site[i,t]] #+ betas[4]*hacked[i] # first year to breeder
logit(psiAB[i,t]) <- eta[5, site[i,t],t] +
lmus[5, site[i,t]] + betas[5]*hacked[i] # nonbreeder to breeder
logit(psiBA[i,t]) <-
lmus[6, site[i,t]] #+ betas[6]*hacked[i] # breeder to nonbreeder
#Re-encounter
logit(pA[i,t]) <- eta[7, site[i,t],t] +
lmus[7, site[i,t]] + betas[7]*hacked[i] +
deltas[5]*effort2[t, site[i,t]] + deltas[6]*effort2[t, site[i,t]]^2 # resight of nonbreeders
logit(pB[i,t]) <- eta[8, site[i,t],t] +
lmus[8, site[i,t]] + #betas[8]*hacked[i] +
deltas[7]*effort2[t, site[i,t]] + deltas[8]*effort2[t, site[i,t]]^2 # resight of breeders
}#t
}#i
# Future projections of survival, recruitment, and detection
cap <- 0.97
for(sc in 1:SC){
for (t in nyr:(nyr+K)){
for (s in 1:nsite){
#Survival
mn.phiFY1[sc, t, s] <- ilogit( eta[1, s, t] +
lmus[1, s] ) + surv.diff[sc, 1, s]
mn.phiA1[sc, t, s] <- ilogit( eta[2, s, t] +
lmus[2, s] ) + surv.diff[sc, 2, s]
mn.phiB1[sc, t, s] <- ilogit( eta[3, s, t] +
lmus[3, s] ) + surv.diff[sc, 3, s]
# enforce the cap
mn.phiFY[sc, t, s] <- step(cap - mn.phiFY1[sc, t, s]) * mn.phiFY1[sc, t, s] +
(1 - step(cap - mn.phiFY1[sc, t, s])) * cap
mn.phiA[sc, t, s] <- step(cap - mn.phiA1[sc, t, s]) * mn.phiA1[sc, t, s] +
(1 - step(cap - mn.phiA1[sc, t, s])) * cap
mn.phiB[sc, t, s] <- step(cap - mn.phiB1[sc, t, s]) * mn.phiB1[sc, t, s] +
(1 - step(cap - mn.phiB1[sc, t, s])) * cap
#Recruitment
logit( mn.psiFYB[sc, t, s] ) <- eta[4, s, t] +
lmus[4, s] # first year to breeder
logit( mn.psiAB[sc, t, s] ) <- eta[5, s, t] +
lmus[5, s] # nonbreeder to breeder
logit( mn.psiBA[sc, t, s] ) <-
lmus[6, s] # breeder to nonbreeder
#Re-encounter
logit( mn.pA[sc, t, s] ) <- eta[7, s, t] +
lmus[7, s] # resight of nonbreeders
logit( mn.pB[sc, t, s] ) <- eta[8, s, t] +
lmus[8, s] # resight of breeders
} # s
} # t
} # sc
# Define state-transition and observation matrices
for (i in 1:nind){
# Define probabilities of state S(t+1) given S(t)
for (t in first[i]:(nyr-1)){
ps[1,i,t,1] <- 0
ps[1,i,t,2] <- phiFY[i,t] * (1-psiFYB[i,t])
ps[1,i,t,3] <- phiFY[i,t] * psiFYB[i,t]
ps[1,i,t,4] <- (1-phiFY[i,t])
ps[2,i,t,1] <- 0
ps[2,i,t,2] <- phiA[i,t] * (1-psiAB[i,t])
ps[2,i,t,3] <- phiA[i,t] * psiAB[i,t]
ps[2,i,t,4] <- (1-phiA[i,t])
ps[3,i,t,1] <- 0
ps[3,i,t,2] <- phiB[i,t] * psiBA[i,t]
ps[3,i,t,3] <- phiB[i,t] * (1-psiBA[i,t])
ps[3,i,t,4] <- (1-phiB[i,t])
ps[4,i,t,1] <- 0
ps[4,i,t,2] <- 0
ps[4,i,t,3] <- 0
ps[4,i,t,4] <- 1
# Define probabilities of O(t) given S(t)
po[1,i,t,1] <- 1
po[1,i,t,2] <- 0
po[1,i,t,3] <- 0
po[1,i,t,4] <- 0
po[2,i,t,1] <- 0
po[2,i,t,2] <- pA[i,t]
po[2,i,t,3] <- 0
po[2,i,t,4] <- (1-pA[i,t])
po[3,i,t,1] <- 0
po[3,i,t,2] <- 0
po[3,i,t,3] <- pB[i,t]
po[3,i,t,4] <- (1-pB[i,t])
po[4,i,t,1] <- 0
po[4,i,t,2] <- 0
po[4,i,t,3] <- 0
po[4,i,t,4] <- 1
} #t
} #i
# Likelihood
for (i in 1:nind){
# Define latent state at first capture
z[i,first[i]] <- y[i,first[i]]
for (t in (first[i]+1):nyr){
# State process: draw S(t) given S(t-1)
z[i,t] ~ dcat(ps[z[i,t-1], i, t-1, 1:4])
# Observation process: draw O(t) given S(t)
y[i,t] ~ dcat(po[z[i,t], i, t-1, 1:4])
} #t
} #i
} )
#**********************
#* Function to run model in NIMBLE
#**********************
run_pva <- function(info, datl, constl, code){
library('nimble')
library('coda')
# function for multivariate normal
uppertri_mult_diag <- nimbleFunction(
run = function(mat = double(2), vec = double(1)) {
returnType(double(2))
p <- length(vec)
out <- matrix(nrow = p, ncol = p, init = FALSE)
for(i in 1:p){
out[ ,i] <- mat[ ,i] * vec[i]
}
return(out)
})
assign('uppertri_mult_diag', uppertri_mult_diag, envir = .GlobalEnv)
params <- c(
# pop growth
"lambda",
# prod
"lmu.prod", "gamma", "rr", "mn.prod",
# survival
"mus", "lmus", "betas", "deltas",
# abundance
"NB", "NF", "NFY", "N", "NAD", "Ntot",
"r",
# error terms
"eta", "sds2", "Ustar2", "U2", "R2",
# yearly summaries
'mn.phiFY', 'mn.phiA', 'mn.phiB',
'mn.psiFYB', 'mn.psiAB', 'mn.psiBA',
'mn.pA', 'mn.pB',
# pva
"extinct", "extinctAD", "extinctB"
)
#n.chains=1; n.thin=1; n.iter=50; n.burnin=25
n.chains=1; n.thin=400; n.iter=1000000; n.burnin=600000
mod <- nimbleModel(code,
constants = constl,
data = datl,
inits = info$inits,
buildDerivs = FALSE,
calculate=T
)
cmod <- compileNimble(mod, showCompilerOutput = TRUE)
confhmc <- configureMCMC(mod)
confhmc$setMonitors(params)
hmc <- buildMCMC(confhmc)
chmc <- compileNimble(hmc, project = mod,
resetFunctions = TRUE,
showCompilerOutput = TRUE)
post <- runMCMC(chmc,
niter = n.iter,
nburnin = n.burnin,
nchains = n.chains,
thin = n.thin,
samplesAsCodaMCMC = T,
setSeed = info$seed)
return(post)
} # run_pva function end
#*****************
#* Run chains in parallel
#*****************
this_cluster <- makeCluster(cpus)
post <- parLapply(cl = this_cluster,
X = par_info_pva[1:cpus],
fun = run_pva,
datl = datl,
constl = constl,
code = mycode)
stopCluster(this_cluster)
save(post, mycode, seeds, cpus,
file="/bsuscratch/brianrolek/riha_ipm/outputs/pva_longrun.rdata")
# save(post, mycode, seeds, cpus,
# file="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//outputs//pva_survival_longrun.rdata")2. Figures and tables
2.1 Plot model estimates
library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
library ("tidybayes")
library ('bayestestR')
library ("ggpubr")
library("viridis")
library ("HDInterval")
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_longrun.rdata")
load("data/data.rdata")
out <- lapply(post, as.mcmc)
outp <- MCMCpstr(out, type="chains")
#********************
#* Calculate summary stats
#********************
#* Abundance at both sites
Nall <- apply(outp$Ntot, c(1,3) , sum)
l.all.post <- array(NA, dim=c(constl$nyr-1, ncol(Nall)) )
for (t in 2:13){
l.all.post[t-1, ] <- Nall[t,] / Nall[t-1,]
}
l.mn <- apply(l.all.post, 2, mean)
# Average annual pop growth
median(l.mn)## [1] 1.020389HDInterval::hdi(l.mn)## lower upper
## 0.9811828 1.0491875
## attr(,"credMass")
## [1] 0.95# 2011 to 2023 pop growth, growth between years
l.all.post2 <- Nall[13,] / Nall[1,]
median(l.all.post2)## [1] 1.246667HDInterval::hdi(l.all.post2)## lower upper
## 0.7111111 1.6995885
## attr(,"credMass")
## [1] 0.95# Population growth rates
mn.lambda <- apply(outp$lambda, c(2,3), mean)
(md.lambda <- apply(mn.lambda, 1, median))## [1] 1.006388 1.298667(hdi95.lambda <- apply(mn.lambda, 1, HDInterval::hdi))## [,1] [,2]
## lower 0.9669408 1.206648
## upper 1.0360502 1.420549# Overall averages
mn.mus <- apply(outp$mus, c(1,3), mean)
(md.mus <- apply(mn.mus, 1, median)) |> round(2)## [1] 0.31 0.94 0.89 0.12 0.43 0.07 0.21 0.77(hdi.mus <- apply(mn.mus, 1, HDInterval::hdi)) |> round(2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.17 0.82 0.82 0.02 0.28 0.03 0.04 0.53
## upper 0.47 1.00 0.96 0.31 0.55 0.12 0.43 0.96# Compare between sites
df.sites <- data.frame(mus=1:8, pd=rep(NA,8), greater=NA)
for (i in 1:8){
first <- median(outp$mus[i, 1, ]) > median(outp$mus[i, 2, ])
df.sites$greater[i] <- ifelse(first, "LH", "PC")
if (first){
diff <- outp$mus[i, 1, ]-outp$mus[i, 2, ]
df.sites$pd[i] <- mean(diff>0) |> round(2)
} else{
diff <- outp$mus[i, 2, ]-outp$mus[i, 1, ]
df.sites$pd[i] <- mean(diff>0) |> round(2)
}
}
df.sites(md.mus <- apply(outp$mus, c(1,2), median)) |> round(2)## [,1] [,2]
## [1,] 0.35 0.26
## [2,] 0.97 0.92
## [3,] 0.91 0.87
## [4,] 0.14 0.07
## [5,] 0.12 0.73
## [6,] 0.12 0.03
## [7,] 0.06 0.33
## [8,] 0.68 0.88(hdi.mus <- apply(outp$mus, c(1,2), HDInterval::hdi)) |> round(2)## , , 1
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.17 0.89 0.84 0.01 0.05 0.05 0.00 0.27
## upper 0.57 1.00 0.98 0.38 0.19 0.20 0.25 0.97
##
## , , 2
##
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## lower 0.08 0.68 0.74 0.00 0.45 0.00 0.04 0.57
## upper 0.48 1.00 0.97 0.38 0.95 0.08 0.73 1.00# productivity
(outp$lmu.prod[1,]-outp$lmu.prod[2,]) |> pd()## Probability of Direction: 1.00# Compare survival, recruitment, and detection rates
# among life stages within each site
df.comp <- data.frame(param1 = c('phiA', 'phiA', 'phiB', 'psiA', 'pB'),
param2 = c('phiB', 'phiFY', 'phiFY', 'psiFY', 'pA'),
pd_LH = NA, pd_PC = NA)
## Los Haitises, survival
df.comp$pd_LH[1] <- (outp$lmus[2,1,]-outp$lmus[3,1,]) |> pd() |> round(2)
df.comp$pd_LH[2] <- (outp$lmus[2,1,]-outp$lmus[1,1,]) |> pd() |> round(2)
df.comp$pd_LH[3] <- (outp$lmus[3,1,]-outp$lmus[1,1,]) |> pd() |> round(2)
## Los Haitises, recruitment rates
df.comp$pd_LH[4] <- (outp$lmus[5,1,]-outp$lmus[4,1,]) |> pd() |> round(2)
## Los Haitises, resight probability
df.comp$pd_LH[5] <- (outp$lmus[8,1,]-outp$lmus[7,1,]) |> pd() |> round(2)
## Punta Canas, survival
df.comp$pd_PC[1] <- (outp$lmus[2,2,]-outp$lmus[3,2,]) |> pd() |> round(2)
df.comp$pd_PC[2] <- (outp$lmus[2,2,]-outp$lmus[1,2,]) |> pd() |> round(2)
df.comp$pd_PC[3] <- (outp$lmus[3,2,]-outp$lmus[1,2,]) |> pd() |> round(2)
## Punta Cana, recruitment rates
df.comp$pd_PC[4] <- (outp$lmus[5,2,]-outp$lmus[4,2,]) |> pd() |> round(2)
## Punta Cana, resight probability
df.comp$pd_PC[5] <- (outp$lmus[8,2,]-outp$lmus[7,2,]) |> pd() |> round(2)
df.comp #*******************
#* Plot IPM results
#*******************
#Plot population size
lNall <- melt(Nall)
colnames(lNall) <- c("Time", "Iter", "Abundance")
lNall$Site.ind <- 3; lNall$Site <- "Both"
lNall <- lNall[,c(1,4,2,3,5)]
Ntot <- outp$Ntot[1:13,,] # thin by 10 to speed up plotting
lN <- melt(Ntot)
colnames(lN) <- c("Time", "Site.ind", "Iter", "Abundance")
lN$Site <- ifelse(lN$Site.ind==1, "Los Haitises", "Punta Cana")
lN <- rbind(lN, lNall)
lN$Year <- lN$Time + 2010
med_hdis <- function(x, cm){
df <- data.frame(y=median(x),
ymin=HDInterval::hdi(x, credMass=cm)[1],
ymax=HDInterval::hdi(x, credMass=cm)[2] )
}
counts.tot <- datl$countsAdults+datl$countsFY+constl$hacked.counts[,-3]
df.counts <- melt(counts.tot)
colnames(df.counts) <- c("Year", "Site_ab", "Count")
df.counts$Site <- ifelse(df.counts$Site_ab=="LHNP", "Los Haitises", "Punta Cana")
df.counts.all <- data.frame(Year=2011:2023, Site_ab="Both", Count=counts.tot[,1] + counts.tot[,2], Site="Both" )
df.counts <- rbind(df.counts, df.counts.all)
p1 <- ggplot() + theme_minimal(base_size=14) +
geom_line(data=lN, aes(x=Year, y=Abundance, group=Iter),
color="gray30", linewidth=0.1, alpha=0.01 ) +
stat_summary(data=lN, aes(x=Year, y=Abundance),
geom="errorbar", width=0,
fun.data =med_hdis, fun.args=list(cm=0.95) ) +
stat_summary(data=lN, aes(x=Year, y=Abundance),
geom="errorbar", width=0, linewidth=1,
fun.data =med_hdis, fun.args=list(cm=0.80) ) +
stat_summary(data=lN, aes(x=Year, y=Abundance),
geom="line", linewidth=1,
fun.data =function(x){data.frame(y=median(x))} ) +
geom_line(data=df.counts, aes(x=Year, y=Count),
col="black", linewidth=1, linetype=2 ) +
facet_wrap("Site", scales="free_y") +
xlab("Year") + ylab("Number")
p1
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Abundance and counts.tiff",
# plot=p1,
# device="tiff",
# width=9, height=4, dpi=300)# Plot life stage structure
mdFY <- apply(outp$NFY[1:13,,], c(1,2), median)
mdB <- apply(outp$NB[1:13,,], c(1,2), median)
mdF <- apply(outp$NF[1:13,,], c(1,2), median)
lFY <- melt(mdFY)
lB <- melt(mdB)
lF <- melt(mdF)
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
ldat <- rbind(lFY, lB, lF)
colnames(ldat)[1:3] <- c("Year", "Sitenum", "Number")
ldat$Site <- ifelse(ldat$Sitenum==1, "Los Haitises", "Punta Cana")
ldat$year <- ldat$Year+2010
# Use median number of females in each stage
# to plot an approximate population structure
p3 <- ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=year)) +
geom_bar(position="fill", stat="identity") +
ylab("Proportion of population") + xlab("Year") +
facet_wrap("Site") +
theme_minimal(base_size=14)+ theme(legend.position="top") +
scale_fill_viridis_d(option = "mako")
p4 <- ggplot(ldat, aes(fill=Stage, y=as.numeric(Number), x=year)) +
geom_bar(position="stack", stat="identity") +
ylab("Median number of females") + xlab("Year") +
facet_wrap("Site", scales = "free_y") +
theme_minimal(base_size=14 ) + theme(legend.position="top") +
scale_fill_viridis_d(option = "mako")
# plot with uncertainty
lFY <- melt(outp$NFY)
lB <- melt(outp$NB)
lF <- melt(outp$NF)
colnames(lFY) <- colnames(lB) <- colnames(lF) <- c("Time", "Site", "Iteration", "Abundance")
lFY$Year <- lFY$Time + 2010
lB$Year <- lB$Time + 2010
lF$Year <- lF$Time + 2010
lFY$Stage <- "First-year"
lB$Stage <- "Breeder"
lF$Stage <- "Nonbreeder"
lALL <- rbind(lFY, lB, lF)
lALL$Site <- ifelse(lALL$Site==1, "Los Haitises", "Punta Cana")
mdFY <- apply(outp$NFY[1:13,,], c(1,2), median)
mdB <- apply(outp$NB[1:13,,], c(1,2), median)
mdF <- apply(outp$NF[1:13,,], c(1,2), median)
hdiFY <- apply(outp$NFY[1:13,,], c(1,2), HDInterval::hdi)
hdiB <- apply(outp$NB[1:13,,], c(1,2), HDInterval::hdi)
hdiF <- apply(outp$NF[1:13,,], c(1,2), HDInterval::hdi)
medFY <- melt(mdFY)
medB <- melt(mdB)
medF <- melt(mdF)
hdisFY <- melt(hdiFY)
hdisB <- melt(hdiB)
hdisF <- melt(hdiF)
dfstages <- rbind(medFY, medB, medF)
dfhdis <- rbind(hdisFY, hdisB, hdisF)
dfstages$LHDI <- dfhdis[dfhdis$Var1=="lower",]
dfstages$UHDI <- dfhdis[dfhdis$Var1=="upper",]
cols <- mako(3)
med_hdis <- function(x, cm){
df <- data.frame(y=median(x),
ymin=HDInterval::hdi(x, credMass=cm)[1],
ymax=HDInterval::hdi(x, credMass=cm)[2] )
}
p5 <- ggplot() + theme_minimal(base_size=14) +
theme(legend.position="none") +
geom_line(data=lALL, aes(x=Year, y=Abundance,
group=Iteration, color=Stage,
alpha=Stage),
linewidth=0.1) +
stat_summary(data=lALL, aes(x=Year, y=Abundance),
geom="errorbar", width=0,
fun.data =med_hdis, fun.args=list(cm=0.95) ) +
stat_summary(data=lALL, aes(x=Year, y=Abundance),
geom="errorbar", width=0, linewidth=1,
fun.data =med_hdis, fun.args=list(cm=0.80) ) +
stat_summary(data=lALL, aes(x=Year, y=Abundance),
geom="line", linewidth=1,
fun.data =function(x){data.frame(y=median(x))} ) +
scale_color_manual( values = c("First-year"=cols[2], "Breeder"=cols[1], "Nonbreeder"=cols[3]) ) +
scale_alpha_manual( values = c("First-year"=0.015, "Breeder"=0.0075, "Nonbreeder"=0.15) ) +
facet_wrap(Stage~Site, scales="free_y",
shrink=TRUE, ncol=2) +
xlab("Year") + ylab("Number of females") +
theme(
strip.background = element_blank(),
strip.text.x = element_blank()
)
p35 <- ggarrange(p3, p5, ncol=1, nrow=2)
p35
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Stage structure_abundance.tiff",
# plot=p35,
# device="tiff",
# width=6, height=8, dpi=300)# plot survival, recruitment, and detection
mus.mat <- outp$mus
dimnames(mus.mat)[[1]] <- c("FY",
"A\nSurvival", "B",
"FY:B",
"A:B\nRecruitment",
"B:A",
"A", "B\nDetection")
dimnames(mus.mat)[[2]] <- c("Los Haitises", "Punta Cana")
lmus <- melt(mus.mat)
colnames(lmus) <- c("Parameter", "Site", "Iter", "value")
# p <- list()
# p[[1]] <- apply(outp$mn.phiFY, c(2, 3), mean)
# p[[2]] <- apply(outp$mn.phiA, c(2, 3), mean)
# p[[3]] <- apply(outp$mn.phiB, c(2, 3), mean)
# p[[4]] <- apply(outp$mn.psiFYB, c(2, 3), mean)
# p[[5]] <- apply(outp$mn.psiAB, c(2, 3), mean)
# p[[6]] <- apply(outp$mn.psiBA, c(2, 3), mean)
# p[[7]] <- apply(outp$mn.pA, c(2, 3), mean)
# p[[8]] <- apply(outp$mn.pB, c(2, 3), mean)
# pl <- lapply(p, melt)
# for (i in 1:8){
# nms <- c('FY', 'A', 'B', 'FY:B', 'A:B', 'B:A', 'pA', 'pB')
# pl[[i]]$nm <- nms[i]
# }
# pl <- do.call(rbind, pl)
# colnames(pl) <- c('Site', 'Iter', 'value', 'nm')
# pl$nm <- factor(pl$nm, levels=c('FY', 'A', 'B', 'FY:B', 'A:B', 'B:A', 'pA', 'pB'))
p6 <- ggplot(data= lmus, aes(x=value, y=Parameter )) +
stat_halfeye(point_interval=median_hdi, .width = c(.95, .80)) +
coord_flip() +
facet_wrap("Site") +
xlab("Probability") + ylab("Parameter") +
theme_bw(base_size=14) +
theme(panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank())
p6
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Survival Recruitment Detection.tiff",
# plot=p5,
# device="tiff",
# width=8, height=4, dpi=300)
# print estimates
mds <- tapply(lmus$value, list(lmus$Parameter, lmus$Site), median) |> round(2)
hdi_LH <- tapply(lmus$value, list(lmus$Parameter, lmus$Site), HDInterval::hdi)[,1]
hdi_LH2 <- do.call(rbind, hdi_LH) |> round(2)
hdi_PC <- tapply(lmus$value, list(lmus$Parameter, lmus$Site), HDInterval::hdi)[,2]
hdi_PC2 <- do.call(rbind, hdi_PC) |> round(2)
df <- data.frame(param= c("phiFY", "phiA", "phiB", "psiFY:B", "psiA:B", "psiB:A", "pA", "pB"),
md_LH=mds[,1], lhdi_LH=hdi_LH2[,1], uhdi_LH=hdi_LH2[,2],
md_PC=mds[,2], lhdi_PC=hdi_PC2[,1], uhdi_PC=hdi_PC2[,2])
df# Plot predicted survival, recruitment, and detection
# with effects from translocation and hacking
# Birds were only hacked from LHNP to PC here
# so we only predict values for PC
betas.pd <- apply(outp$betas, 1, pd)
ind <- which (betas.pd>0.975) # which are significant
pred.mus <- array(NA, dim=c(dim(outp$mus)))
for (m in ind){
for (tr in 1:2){
pred.mus[m,tr,] <- plogis( outp$lmus[m,2,] + outp$betas[m,]*c(0,1)[tr] )
}}
# treatment, mu, site, iter
mus.md <- apply(pred.mus, c(1,2), median)
mus.HDI80 <- apply(pred.mus, c(1,2), HDInterval::hdi, credMass=0.8)
mus.HDI95 <- apply(pred.mus, c(1,2), HDInterval::hdi)
# tiff(height=4, width=6, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Translocation effects.tiff")
par(mar=c(3,4,4,1), mfrow=c(1,1))
for (tr in 1:2){
if (tr==1){
plot(1:4, c(0,0,1,1),
type="n",
pch=c(16, 17)[tr], cex=2,
ylim=c(0,1),
xlim=c(0.5, 4.5),
cex.axis=1, cex.lab=1.25,
ylab="Probability", xlab="",
main="",
xaxt="n", yaxt="n")
abline(h= seq(0,1,by=0.1), col="gray90")
abline(v= seq(0.5,6,by=1), col="gray90")
points(1:4+c(-0.1,0.1)[tr], mus.md[ind,tr],
pch=c(16, 17)[tr], cex=1.5)
}else{
points(1:4+c(-0.1,0.1)[tr], mus.md[ind,tr],
pch=c(16, 15)[tr], cex=1.5)
}}
axis(1, at=1:4, cex.axis=0.85,
labels=c("First-year\nSurvival",
"Breeder\nSurvival",
"Nonbreeder\n to Breeder\nRecruitment",
"Nonbreeder\nDetection"), las=1, padj=0.5)
axis(2, at=seq(0,1, by=0.1), cex.axis=1,
labels=c(0, NA,NA,NA,NA,0.5, NA,NA,NA,NA, 1))
for (m in 1:length(ind)){
for (tr in 1:2){
lines(rep(c(1:4)[m] + c(-0.1,0.1)[tr],2), mus.HDI95[1:2,ind[m],tr], lwd=2)
lines(rep(c(1:4)[m] + c(-0.1,0.1)[tr],2), mus.HDI80[1:2,ind[m],tr], lwd=4)
}}
legend(x=2.5,y=1.25, pch=c(16,15), pt.cex=1.5, cex=1, xpd=NA, horiz=T,
legend=c("Not translocated", "Translocated" ),
xjust=0.5)
#dev.off()# Plot fecundity in response to nest treatments
f <- exp(outp$lmu.prod)
f.md <- apply(f, 1, median)
f.HDI80 <- apply(f, 1, HDInterval::hdi, credMass=0.8)
f.HDI95 <- apply(f, 1, HDInterval::hdi)
# Calculate treatment effects on fecundity
f.pred <- array(NA, dim=dim(outp$lmu.prod))
for (s in 1:2){
f.pred[s,] <- exp(outp$lmu.prod[s,] + outp$gamma[,1])
} # s
f2.md <- apply(f.pred, 1, median)
f2.HDI80 <- apply(f.pred, 1, HDInterval::hdi, credMass=0.8)
f2.HDI95 <- apply(f.pred, 1, HDInterval::hdi)
# tiff(height=4, width=6, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Fecundity_treatment.tiff")
par(mar=c(3,5,3,1), mfrow=c(1,1))
plot(c(1, 2)-0.1, f.md,
ylim=c(0, 1.55),
xlim=c(0.5, 2.5),
ylab="Productivity", xlab="",
main="",
cex.axis=1.5, cex.lab=1.5,
xaxt="n", yaxt="n", type="n")
points(c(1, 2)-0.1, f.md, pch=16, cex=1.5)
abline(h=seq(0, 1.6, by=0.1), col="gray90")
abline(v=1.5, col="black")
axis(1, at=c(1,2), cex.axis=1,
labels=c("Los Haitises\nNational Park","Punta Cana"),
padj=0.5)
axis(2, at=seq(0, 1.6, by=0.1),
cex.axis=1,
labels=c("0", NA, NA, NA, NA, "0.5", NA, NA, NA, NA, "1.0", NA, NA, NA, NA, "1.5", NA))
lines(c(1,1)-0.1, f.HDI95[,1], lwd=2)
lines(c(2,2)-0.1, f.HDI95[,2], lwd=2)
lines(c(1,1)-0.1, f.HDI80[,1], lwd=4)
lines(c(2,2)-0.1, f.HDI80[,2], lwd=4)
points(c(1.1, 2.1), f2.md,
pch=18, cex=2.5)
lines(c(1,1) +0.1, f2.HDI95[,1], lwd=3)
lines(c(2,2) +0.1, f2.HDI95[,2], lwd=3)
lines(c(1,1) +0.1, f2.HDI80[,1], lwd=6)
lines(c(2,2) +0.1, f2.HDI80[,2], lwd=6)
legend(x=1.5,y=1.9, pch=c(16,18), pt.cex=c(1.5,2.5), cex=1,
legend=c("Untreated", "Treated" ),
horiz=TRUE, xjust=0.5, xpd=NA)
#dev.off()
# Calculate probability of direction
f.diff <- f[1,]-f[2,]
pd(f.diff)## Probability of Direction: 1.00f.md ## lmu.prod[1] lmu.prod[2]
## 0.3743717 0.2265712f.HDI95 ## lmu.prod[1] lmu.prod[2]
## lower 0.2952690 0.1640874
## upper 0.4527469 0.2986755f.diff2 <- f.pred[1,]-f.pred[2,]
pd(f.diff2)## Probability of Direction: 1.00f2.md ## [1] 1.2559463 0.7601034f2.HDI95## [,1] [,2]
## lower 0.9905714 0.550482
## upper 1.5188800 1.002000median(f.pred[1,]/f[1,])## [1] 3.354811median(f.pred[2,]/f[2,])## [1] 3.354811#*******************
# Correlations between demographic rates
#*******************
pdR2 <- apply(outp$R2, c(1,2), pd) |> round(2)
max(pdR2[pdR2<1])## [1] 0.82apply(outp$R2, c(1,2), median) |> round(2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1.00 -0.02 0.14 0.11 -0.07 0.00 0.08 0.07 0.05
## [2,] -0.02 1.00 -0.02 -0.03 0.02 0.00 -0.07 -0.02 0.00
## [3,] 0.14 -0.02 1.00 -0.05 -0.11 0.01 0.20 -0.06 -0.04
## [4,] 0.11 -0.03 -0.05 1.00 -0.11 0.00 0.23 -0.22 -0.17
## [5,] -0.07 0.02 -0.11 -0.11 1.00 -0.01 -0.08 -0.02 0.08
## [6,] 0.00 0.00 0.01 0.00 -0.01 1.00 0.01 0.00 0.00
## [7,] 0.08 -0.07 0.20 0.23 -0.08 0.01 1.00 -0.22 -0.21
## [8,] 0.07 -0.02 -0.06 -0.22 -0.02 0.00 -0.22 1.00 0.04
## [9,] 0.05 0.00 -0.04 -0.17 0.08 0.00 -0.21 0.04 1.00apply(outp$R2, c(1,2), HDInterval::hdi)[1,,] |> round(2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1.00 -0.64 -0.45 -0.50 -0.64 -0.55 -0.46 -0.45 -0.52
## [2,] -0.64 1.00 -0.61 -0.62 -0.57 -0.57 -0.63 -0.59 -0.59
## [3,] -0.45 -0.61 1.00 -0.63 -0.66 -0.62 -0.35 -0.57 -0.54
## [4,] -0.50 -0.62 -0.63 1.00 -0.68 -0.57 -0.34 -0.71 -0.71
## [5,] -0.64 -0.57 -0.66 -0.68 1.00 -0.62 -0.61 -0.57 -0.49
## [6,] -0.55 -0.57 -0.62 -0.57 -0.62 1.00 -0.60 -0.62 -0.63
## [7,] -0.46 -0.63 -0.35 -0.34 -0.61 -0.60 1.00 -0.63 -0.65
## [8,] -0.45 -0.59 -0.57 -0.71 -0.57 -0.62 -0.63 1.00 -0.45
## [9,] -0.52 -0.59 -0.54 -0.71 -0.49 -0.63 -0.65 -0.45 1.00apply(outp$R2, c(1,2), HDInterval::hdi)[2,,] |> round(2)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 1.00 0.56 0.67 0.62 0.51 0.62 0.62 0.59 0.60
## [2,] 0.56 1.00 0.55 0.58 0.59 0.60 0.53 0.58 0.59
## [3,] 0.67 0.55 1.00 0.53 0.48 0.58 0.69 0.45 0.44
## [4,] 0.62 0.58 0.53 1.00 0.44 0.60 0.71 0.31 0.40
## [5,] 0.51 0.59 0.48 0.44 1.00 0.56 0.48 0.52 0.61
## [6,] 0.62 0.60 0.58 0.60 0.56 1.00 0.57 0.59 0.59
## [7,] 0.62 0.53 0.69 0.71 0.48 0.57 1.00 0.23 0.33
## [8,] 0.59 0.58 0.45 0.31 0.52 0.59 0.23 1.00 0.50
## [9,] 0.60 0.59 0.44 0.40 0.61 0.59 0.33 0.50 1.00#*******************
# Correlations between demographics and population growth rates
#*******************
lam.m <- apply(outp$lambda[1:12,,], c(1,2), median)
lam.hdi <- apply(outp$lambda[1:12,,], c(1,2), HDInterval::hdi)
par(mfrow=c(1,2))
plot(2012:2023, lam.m[,1], type="b", pch=1,
ylab="Population growth rate", xlab="Year",
ylim=c(min(lam.hdi[,,1]), max(lam.hdi[,,1])),
main="Los Haitises")
abline(h=1, lty=2)
segments(x0=2012:2023, x1=2012:2023,
y0 = lam.hdi[1,,1], y1= lam.hdi[2,,1])
plot(2012:2023, lam.m[,2], type="b", pch=1, lty=1,
ylab="Population growth rate", xlab="Year",
ylim=c(min(lam.hdi[,,2]), max(lam.hdi[,,2])),
main="Punta Cana")
abline(h=1, lty=2)
segments(x0=2012:2023, x1=2012:2023,
y0 = lam.hdi[1,,2], y1= lam.hdi[2,,2])
# create a function to plot correlations
# between demographics and population growth rates
plot.cor <- function (lam.post, x.post, x.lab, ind.x=1:12, site=1, yaxt="b"){
# calculate the correlation coefficient
# over each iteration to propagate error
lambda <- lam.post[1:12, site,]
x <- x.post[ind.x, site, ]
df <- data.frame(ats1=c(0.8, 0.9, 1.0, 1.1, 1.2, 1.3),
ats2=c(1.0, 1.5, 2.0, 2.5, 3.0, 3.5),
labs1=c("0.8", NA, "1.0", NA, "1.2", NA),
labs2=c("1.0", NA, "2.0", NA, "3.0", NA)
)
if(site==1){ df <- df[, c(1,3)]} else{
df <- df[, c(2,4)]
}
lwd <- 1
cor.post <- array(NA, dim=dim(lambda)[-1])
for (i in 1:dim(lambda)[[2]]){
cor.post[i] <- cor(lambda[,i], x[,i])
}
cor.df <- data.frame(median= median(cor.post) |> round(digits=2),
ldi=HDInterval::hdi(cor.post)[1] |> round(digits=2),
hdi=HDInterval::hdi(cor.post)[2] |> round(digits=2),
pd = pd(cor.post) |> round(digits=2))
lam.m <- apply(lambda, 1, median)
lam.hdi <- apply(lambda, 1, HDInterval::hdi)
x.m <- apply(x, 1, median)
x.hdi <- apply(x, 1, HDInterval::hdi)
x.lims <- c(min(x.hdi[,]), max(x.hdi[,]))
y.lims <- c(min(lam.hdi[,]), max(lam.hdi[,]))
if(yaxt=="n"){
plot(x.m, lam.m[1:12],
xlim= x.lims,
ylim= y.lims,
type="n",
ylab="", xlab=x.lab, cex.lab=1.5,
main=c("", "")[site],
yaxt=yaxt, xaxt="n")
axis(2, at=df[,1], labels=c(rep(NA, 5)))
} else {
plot(x.m, lam.m[1:12],
xlim= x.lims,
ylim= y.lims,
type="n",
ylab="", xlab=x.lab, cex.lab=1.5,
main=c("", "")[site],
yaxt="n", xaxt="n")
axis(2, at=df[,1], labels=df[,2], las=1, cex.axis=1.5)
}
axis(1, cex.axis=1.5)
points(x.m, lam.m, pch=1)
segments(x0=x.hdi[1,], x1=x.hdi[2,],
y0 = lam.m, y1= lam.m, lwd=lwd)
segments(x0=x.m, x1=x.m,
y0 = lam.hdi[1,], y1= lam.hdi[2,], lwd=lwd)
xs <- c(x.lims[1], x.lims[1]+(x.lims[2]-x.lims[1]*1.5))
ys <- c( (y.lims[2]-y.lims[1])+y.lims[1], (y.lims[2]-y.lims[1])+y.lims[1],
(y.lims[2]-y.lims[1])*0.6+y.lims[1], (y.lims[2]-y.lims[1])*0.6+y.lims[1])
polygon(x=c(xs, rev(xs)), y=ys, col=alpha("white", 0.7), border=NA )
text(x = x.lims[1], y = (y.lims[2]-y.lims[1])*0.9+y.lims[1],
paste("r = ", cor.df$median, " (",cor.df$ldi,", ", cor.df$hdi, ")", sep=""),
pos = 4, font = 3, cex = 1.5)
text(x = x.lims[1], y = (y.lims[2]-y.lims[1])*0.7+y.lims[1], paste("P(r>0) = ", cor.df$pd, sep=""),
pos = 4, font = 3, cex = 1.5)
}
# Plor correlations between population grrowth rates
# and demographics.
# "r" is a correlation coefficient and represents
# the magnitude of the correlation
# P(r>0) is the probability of direction (similar to p-values)
# that is, the probability that an effect exists
# tiff(height=8, width=8, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//PopGrowthand Demographics1.tiff")
# par(mfrow=c(4,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
par(mfrow=c(2,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
plot.cor(outp$lambda, outp$mn.prod, x.lab="Fecundity", ind.x=2:13)
plot.cor(outp$lambda, outp$mn.phiFY, x.lab="First-year Survival", yaxt="n")
plot.cor(outp$lambda, outp$mn.phiA, x.lab="Nonbreeder Survival", yaxt="n")
plot.cor(outp$lambda, outp$mn.phiB, x.lab="Breeder Survival")
plot.cor(outp$lambda, outp$mn.psiFYB, x.lab="First-year to Breeder", yaxt="n")
plot.cor(outp$lambda, outp$mn.psiAB, x.lab="Nonbreeder to Breeder", yaxt="n")
mtext("Population growth rate", side=2, outer=TRUE, line=2.5, cex=2)
mtext("Los Haitises", side=3, outer=TRUE, cex=2)
#dev.off()
# tiff(height=4, width=8, units="in", res=300,
# filename= "C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//PopGrowthand Demographics2.tiff")
par(mfrow=c(2,3), mar=c(4,1,1,1), oma=c(0,5,3,0))
plot.cor(outp$lambda, outp$mn.prod, x.lab="Fecundity", ind.x=2:13, site=2)
plot.cor(outp$lambda, outp$mn.phiFY, x.lab="First-year Survival", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.phiA, x.lab="Nonbreeder Survival", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.phiB, x.lab="Breeder Survival", site=2)
plot.cor(outp$lambda, outp$mn.psiFYB, x.lab="First-year to Breeder", yaxt="n", site=2)
plot.cor(outp$lambda, outp$mn.psiAB, x.lab="Nonbreeder to Breeder", yaxt="n", site=2)
mtext("Population growth rate", side=2, outer=TRUE, line=2.5, cex=2)
mtext("Punta Cana", side=3, outer=TRUE, cex=2)
#dev.off()#*******************
#* Plot PVA results
#*******************
library ('viridis')
library ("tidybayes")
library ('coda')
library ('MCMCvis')
library ('reshape2')
library('ggplot2')
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\pva_longrun.rdata")
out <- lapply(post, as.mcmc)
outp <- MCMCpstr(out, type="chains")
scen <- data.frame('Scenarios' = 1:36,
'Number translocated' = rep(c(0,5,10), each=12),
'Number of nests treated' = rep( rep(c(0, 15, 30, 45), each=3), 3),
'Mortality reduction' = rep(c('100% Mortality', '80% Mortality', '60% Mortality'), 12) )
#*********************
#* Line plot of extinction probability over time
#*********************
pe <- apply(outp$extinct, c(1,2,3), mean)
lpe <- melt(pe)
colnames(lpe) <- c('Scenarios', 'Time', 'Site2', 'Extinct')
pe.df <- merge(lpe, scen, by='Scenarios', all.x=TRUE)
pe.df$Year <- pe.df$Time+2023
pe.df$Site <- ifelse(pe.df$Site2==1, "Los Haitises", "Punta Cana")
pe.df[(pe.df$Site2==2 & pe.df$Number.translocated>0),]$Extinct <- NA
pe.df$Mortality.reduction <- factor(pe.df$Mortality.reduction,
levels= c('100% Mortality', '80% Mortality', '60% Mortality'))
p6 <- ggplot(pe.df, aes(x=Year, y=Extinct, group=Scenarios,
color = factor(Number.translocated),
linetype = factor(Number.of.nests.treated))) +
geom_line( linewidth = 1) +
facet_grid(rows= vars(Site), cols= vars(Mortality.reduction),
scales='free_y') +
#facet_wrap(~ Site + Mortality.reduction) +
scale_color_viridis_d(option="viridis", direction=1,
name="Number translocated") +
scale_linetype(name="Number of nests treated") +
scale_x_continuous( name="Future year",
breaks=c(2020, 2030, 2040, 2050, 2060, 2070),
labels=c('', '2030', '', '2050', '', '2070'),
minor_breaks=NULL) +
theme_minimal() +
ylab("Extinction probability")
p6
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Extinction_lines.tiff",
# plot=p6,
# device="tiff",
# width=6, height=4, dpi=300)
#********************
#* PVA heatmap of extinction probability projected 50 years
#********************
pe.df2 <- pe.df[pe.df$Time==50, ]
p7 <- ggplot(pe.df2, aes(x = as.factor(Number.translocated),
y = as.factor(Number.of.nests.treated))) +
geom_tile( aes(fill = Extinct ) ) +
geom_text(aes(label=round(Extinct,2)), color="gray70") +
scale_fill_viridis_c(option="plasma", direction=1,
name = "Extinction\nprobability after\n50 years") +
facet_wrap(~ Site + factor(Mortality.reduction,
levels=c("100% Mortality", "80% Mortality", "60% Mortality"))) +
scale_y_discrete(breaks=c(0,15,30,45),
name="Number of nests treated") +
scale_x_discrete(name="Number translocated") +
theme_minimal( ) + theme(strip.background = element_rect(size = 0.5))
p7 
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Extinction_50yrHeatmap.tiff",
# plot=p7,
# device="tiff",
# width=6, height=4, dpi=300)
pe.df3 <- pe.df2[!is.na(pe.df2$Extinct),]
pe.df3 <- pe.df3[rev(order(pe.df3$Site2, pe.df3$Extinct)),]
pe.df3#***************************
# Plot total Abundance
lnt1 <- outp$Ntot[,1:13, , ] |> melt()
colnames(lnt1) <- c("Scenarios", "Time", "Site2", "Iter", "Abund")
lnt1$Year <- lnt1$Time + 2010
lnt1$Site <- factor(ifelse(lnt1$Site2==1, 'Los Haitises', 'Punta Cana'))
nt.df <- merge(lnt1, scen, by='Scenarios', all.x=TRUE)
nt.df$Mortality.reduction <- factor(nt.df$Mortality.reduction,
levels=c('100% Mortality','80% Mortality', '60% Mortality' ))
nt.df[(nt.df$Site2==2 & nt.df$Number.translocated>0),]$Abund <- NA
nt.df <- nt.df[nt.df$Scenarios %in% c(1,2,3), ]
# calculate future means
lnt2 <- apply(outp$Ntot[,13:63, , ], c(1, 2, 3), median) |> melt()
colnames(lnt2) <- c("Scenarios", "Time", "Site2", "Abund")
lnt2$Year <- lnt2$Time + 2022
lnt2$Site <- factor(ifelse(lnt2$Site2==1, 'Los Haitises', 'Punta Cana'))
nt.df2 <- merge(lnt2, scen, by='Scenarios', all.x=TRUE)
nt.df2$Mortality.reduction <- factor(nt.df2$Mortality.reduction,
levels=c('100% Mortality','80% Mortality', '60% Mortality' ))
nt.df2[(nt.df2$Site2==2 & nt.df2$Number.translocated>0),]$Abund <- NA
p8 <- ggplot() +
stat_lineribbon(data= nt.df,
aes(x=Year, y=Abund, group=Scenarios),
fill='gray40',
linewidth=0.5,
point_interval = 'median_hdci', .width = 0.95,
na.rm=TRUE) +
geom_line(data= nt.df2,
aes(x=Year, y=Abund, group=Scenarios,
color = factor(Number.translocated),
linetype = factor(Number.of.nests.treated) ),
linewidth=0.5) +
facet_grid(rows = vars(Site), cols = vars(Mortality.reduction),
scales ='free_y') +
ylab('Number of females') +
scale_color_viridis_d(option="viridis", direction=1,
name="Number translocated") +
scale_linetype(name="Number of nests treated") +
theme_minimal()
p8
# ggsave(filename="C://Users//rolek.brian//OneDrive - The Peregrine Fund//Documents//Projects//Ridgways IPM//figs//Abundance_projections.tiff",
# plot=p9,
# device="tiff",
# width=8, height=4, dpi=300)
# Time to extinction
# conditional on extinction threshold (D=X)D <- 1
abund <- outp$Ntot[,14:(13+50),,]
dimnames(abund) <- list('Scenarios'=1:36,
'Time'= 2024:(2024+49),
'Site'=c("Los Haitises", "Punta Cana"),
'Iter'=1:10000)
qextinct <- abund == D
ext <- qextinct[,50,,] <= D
lte <- melt( qextinct )
colnames(lte)[5] <- "Extinct"
le <- melt( ext )
colnames(le)[4] <- "EverExtinct"
lte2 <- merge(lte, le, by=c("Scenarios", "Site", "Iter") )
lte3<- lte2[lte2$EverExtinct==TRUE, ]
lte4 <- lte3[lte3$Extinct==TRUE, ]
lte5 <- tapply(lte4$Time-2023, list(lte4$Scenarios, lte4$Site, lte4$Iter), min)
lte6 <- melt(lte5)
lte7 <- lte6[!is.na(lte6$value), ]
colnames(lte7) <- c('Scenarios', 'Site', 'Iter', 'TTE')
te <- table(lte7$TTE, lte7$Scenarios, lte7$Site) |> melt()
colnames(te) <- c('TE', 'Scen', 'Site', 'Count')
te$Scenarios <- factor(paste0("Scenario ", te$Scen),
levels=paste0('Scenario ', 1:36))
teLH <- te[te$Site=="Los Haitises",]
tePC <- te[te$Site=="Punta Cana",]
p9 <- ggplot() +
geom_col(data = teLH, aes(x= TE, y=Count)) +
facet_wrap(vars(Scenarios),
nrow=12, ncol=3, scales="free_y", shrink=FALSE)
p10 <- ggplot() +
geom_col(data = tePC, aes(x= TE, y=Count)) +
facet_wrap(vars(Scenarios),
nrow=12, ncol=3, scales="free_y", shrink=FALSE)
p9
p10
2.2 Plots of finer population segments and model diagnostics
load("C:\\Users\\rolek.brian\\OneDrive - The Peregrine Fund\\Documents\\Projects\\Ridgways IPM\\outputs\\ipm_longrun.rdata")
load("data/data.rdata")
library ('MCMCvis')
library ('coda')
library ('ggplot2')
library('reshape2')
library('bayestestR')
out <- lapply(post, as.mcmc)
# Identify chains with NAs that
# failed to initialize
NAlist <- c()
for (i in 1:length(out)){
NAlist[i] <- any (is.na(out[[i]][,1:286]) | out[[i]][,1:286]<0)
}
# Subset chains to those with good initial values
out <- out[!NAlist]
post2 <- post[!NAlist]
outp <- MCMCpstr(out, type="chains")
!NAlist## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE# default settings for plots
plt <- function(object, params,...) {
MCMCplot(object=out,
params=params,
guide_axis=TRUE,
HPD=TRUE, ci=c(80, 95), horiz=FALSE,
#ylim=c(-10,10),
...)
}Plot model estimates of demographic rates. Life Stages are abbreviated as B = breeder, NB = nonbreeder, FY = first year. First-year abundance accounts for translocated birds.
# Abundance of females at Los Haitises
par(mfrow=c(5,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NFY[",1:13, ", 1]"),
main="First-year (FY) minus hacked\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,1],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NF[",1:13, ", 1]"),
main="Adult nonbreeder (NB)\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NB[",1:13, ", 1]"),
main="Adult breeder (B)\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NAD[",1:13, ", 1]"),
main="Adult Breeders and Nonbreeders\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsAdults[,1],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("Ntot[",1:13, ", 1]"),
main="All stages\n Los Haitises",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,1]+datl$countsAdults[,1],
ylab="Counts", xlab="Year", type="b")
# Abundance of females at Punta Cana
par(mfrow=c(5,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NFY[",1:13, ", 2]"),
main="First-year (FY) plus hacked\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,2],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NF[",1:13, ", 2]"),
main="Adult nonbreeder (NB)\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NB[",1:13, ", 2]"),
main="Adult breeder (B)\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot.new()
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("NAD[",1:13, ", 2]"),
main="Adult Breeders and Nonbreeders\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsAdults[,2],
ylab="Counts", xlab="Year", type="b")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("Ntot[",1:13, ", 2]"),
main="All stages\n Punta Cana",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plot(2011:2023, datl$countsFY[,2]+datl$countsAdults[,2],
ylab="Counts", xlab="Year", type="b")
Population dynamics are determined by transitions, These transitions between stages are abbreviated as the starting life stage to the final life stage. For example a first-year recruiting to a breeder would be abbreviated as “FY to B”. I’ll list them here for convenience:
“FY to NB” is first-year to nonbreeder.
“NB to NB” is nonbreeder adult to nonbreeder adult.
“B to NB” is a breeding adult to a nonbreeder adult.
“FY to B” is first-year to breeder.
“NB to B” is nonbreeder adult to breeder adult.
“B to B” is breeder adult to breeder adult.
# Finer population segments
par(mfrow=c(4,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 1, ", ", 1:13, ", 1]"),
main="Los Haitises\nFirst-years fledged",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 2, ", ", 1:13, ", 1]"),
main="Los Haitises\nFY to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 3, ", ", 1:13, ", 1]"),
main="Los Haitises\nNB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 4, ", ", 1:13, ", 1]"),
main="Los Haitises\nB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 5, ", ", 1:13, ", 1]"),
main="Los Haitises\nFY to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 6, ", ", 1:13, ", 1]"),
main="Los Haitises\nNB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 7, ", ", 1:13, ", 1]"),
main="Los Haitises\nB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
par(mfrow=c(4,2))
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 1, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY fledged",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 2, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 3, ", ", 1:13, ", 2]"),
main="Punta Cana\nNB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 4, ", ", 1:13, ", 2]"),
main="Punta Cana\nB to NB",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 5, ", ", 1:13, ", 2]"),
main="Punta Cana\nFY to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 6, ", ", 1:13, ", 2]"),
main="Punta Cana\nNB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
plt(object=out,
exact=TRUE, ISB=FALSE,
params=paste0("N[", 7, ", ", 1:13, ", 2]"),
main="Punta Cana\nB to B",
labels = 2011:2023,
xlab = "Year", ylab= "Abundance")
Other parameter estimates.
# I needed to abbreviate to save plot space
# FY=first-year, NB=Nonbreeder, B=Breeder
par(mfrow=c(1,2))
plt(object=out,
params=paste0("mus[",1:8, ", 1]"),
exact=TRUE, ISB=FALSE,
ylim=c(0,1),
main="Overall means\n Los Haitises",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection")
)
plt(object=out,
params=paste0("mus[",1:8, ", 2]"),
exact=TRUE, ISB=FALSE,
ylim=c(0,1),
main="Overall means\n Punta Cana",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection"))
par(mfrow=c(1,1))
plt(object=out,
params="betas",
main= "Translocation effects",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection"))
# Plot effort effects
plt(object=out,
params=paste0("deltas[",1:4,"]"),
exact=TRUE, ISB=FALSE,
main= "Effort effects",
labels=c("Number of Adults\nEffort", "Number of Adults\nEffort squared",
"Number of FYs\nEffort", "Number of FYs\nEffort squared"),
ylim=c(-0.5, 1))
plt(object=out,
params=paste0("deltas[",5:8,"]"),
exact=TRUE, ISB=FALSE,
main= "Effort effects continued",
labels=c("Nonbreeder detection\nEffort", "Nonbreeder detection\nEffort sq",
"Breeder detection\nEffort", "Breeder detection\nEffort sq"),
ylim=c(-5, 5))
# Fecundity
par(mfrow=c(1,1))
plt(object=out,
params=c("lmu.prod"),
labels= c("Productivity\n(log scale)\nLos Haitises",
"Productivity\n(log scale)\nPunta Cana"))
f <- exp(outp$lmu.prod)
# Is fecundity at LHNP greater than PC
par(mfrow=c(1,1))
fdiff <- f[1,]-f[2,]
hist(fdiff, main="Productivity difference")
abline(v=0, lty=2)
# print probability of direction, similar to frequentist p-value
# so values <=0.025 and >=0.975
# Is the difference in fecundity >0 ?
mean(fdiff>0) ## [1] 0.9988# How many times greater is fecundity
# at treated versus non-treated sites
# median(f.pred[1,]/f[1,])
# median(f.pred[2,]/f[2,])
# gamma = nest treatment effect on fecundity
par(mfrow=c(1,1))
plt(object=out,
params=c("gamma"),
main="Anti-Parasitic Fly\nTreatment Effects", ylim=c(0,3))
par(mfrow=c(1,1))
# sds <- paste0("sds[", 1:9, "]")
# plt(object=out, params=sds,
# exact=TRUE, ISB=FALSE,
# main="Temporal SDs (synchrony among sites)",
# labels=c("FY survival", "NB survival", "B survival",
# "FY to B", "NB to B", "B to NB",
# "NB detection", "B detection",
# "Productivity"))
sds2 <- paste0("sds2[", 1:9, "]")
plt(object=out, params=sds2,
exact=TRUE, ISB=FALSE,
main="Site-temporal SDs",
labels=c("FY survival", "NB survival", "B survival",
"FY to B", "NB to B", "B to NB",
"NB detection", "B detection",
"Productivity"))
# Correlations among vital rates
# Plot is messy with only a few strong correlations
ind <- 1
Rs <- R2s <- c()
for (i in 1:(nrow(outp$R)-1)){
for (j in (i+1):nrow(outp$R)){
Rs[ind] <- paste0("R[",i,", ", j, "]")
R2s[ind] <- paste0("R2[",i,", ", j, "]")
ind <- ind+1
}}
# par(mfrow=c(2,1))
# plt(object=out, params=Rs[1:18], exact=TRUE, ISB=FALSE,
# main="Correlations btw demographic rates\n over time (synchrony)",
# xlab = "Rhos", guide_lines=TRUE)
# plt(object=out, params=Rs[19:36], exact=TRUE, ISB=FALSE,
# main="Correlations btw demographic rates\n over time (synchrony), continued...",
# xlab = "Rhos", guide_lines=TRUE)
par(mfrow=c(2,1))
plt(object=out, params=R2s[1:18], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time and sites",
xlab = "Rhos", guide_lines=TRUE)
plt(object=out, params=R2s[19:36], exact=TRUE, ISB=FALSE,
main="Correlations btw demographic rates\n over time and sites, continued ...",
xlab = "Rhos", guide_lines=TRUE)
# Annual averages for integration into the population model
par(mfrow=c(1,1))
labs <- c(paste0("LH ",2011:2023), paste0("PC ",2011:2023))
plt(object=out, params="mn.phiFY", ylim=c(0,1),
main="First-year survival", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.phiA", ylim=c(0,1),
main="Adult nonbreeder", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.phiB", ylim=c(0,1),
main="Breeder", labels = labs,
xlab = "Year", ylab= "Survival")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiFYB", ylim=c(0,1),
main="First-year to breeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiAB", ylim=c(0,1),
main="Adult nonbreeder to breeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.psiBA", ylim=c(0,1),
main="Adult breeder to nonbreeder", labels = labs,
xlab = "Year", ylab= "Recruitment")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.pA", ylim=c(0,1),
main="Nonbreeder", labels = labs,
xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.pB", ylim=c(0,1),
main="Breeder", labels = labs,
xlab = "Year", ylab= "Detection")
abline(v=13.5, lwd=2)
plt(object=out, params="mn.prod",
main="", labels=labs,
xlab = "Year", ylab= "Productivity")
abline(v=13.5, lwd=2)
2.3 Print parameter estimates
Parameter estimates for input into a population viability analysis.
# Estimates for the survival model
# In this order: FY survival, NB survival, B survival,
# FY to B recruitment, NB to B recruitent, B to NB recruitment,
# Detection NB, Detection B
# Mus are means for
# mus[1, site] , where site=1 is LH and site=2 is PC
# Survival of first years = mus[1,]
# Survival of nonbreeders = mus[2,]
# Survival of breeders = mus[3,]
# Breeding propensity of first years = mus[4,]
# Breeding propensity of nonbreeders = mus[5,]
# Transition from breeder to nonbreeder = mus[6,]
# Detection probability of nonbreeders = mus[7,]
# Detection probability of breeders
# modeled as logit(probability) = lmus[x,site] + betas[x]*translocated + eps[x,t] + eta[x,s,t]
# except breeder to nonbreeder recruitment = lmus[6,site]
MCMCsummary(post2, params = "mus",
digits=3, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")MCMCsummary(post2, params = "betas",
digits=3, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Fecundity
# modeled as log(f) = lmu.f[site] + gamma*treatment + eps[x,t] + eta[x,s,t]
# log scale
MCMCsummary(post2, params = "lmu.prod",
digits=3, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Effort effects on detections and resighting
MCMCsummary(post2, params = "deltas",
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Random effects for temporal synchrony and site x year.
MCMCsummary(post2, params = sds2,
exact=TRUE, ISB=FALSE,
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")# Correlations among demographic rates site x time
MCMCsummary(post2, params = R2s,
exact=TRUE, ISB=FALSE,
digits=2, HPD = T,
hpd_prob = 0.95, pg0= TRUE,
func=median, func_name="median")3. Model diagnostics
3.1 Check Goodness-of-fit
Goodness-of-fit tests provide evidence that statistical distributions adequately describe the data. Here we test fit for brood size and counts. A Bayesian p-value nearest to 0.5 suggests good fitting statistical distributions, while values near 1 or 0 suggest poor fit.
# Goodness of fit check
fit.check <- function(out, ratio=FALSE,
name.rep="f.dmape.rep",
name.obs="f.dmape.obs",
jit=100,
ind=1,
lab=""){
par(mfrow=c(1,1))
# plot mean absolute percentage error
samps <- MCMCpstr(out, "all", type="chains")
rep <- samps[name.rep][[1]][ind,]
obs <- samps[name.obs][[1]][ind,]
mx <- max(c(rep, obs))
mn <- min(c(rep, obs))
plot(jitter(obs, amount=jit),
jitter(rep, amount=jit),
main=paste0("Mean absolute percentage error\n",lab),
ylab="Discrepancy replicate values",
xlab="Discrepancy observed values",
xlim=c(mn,mx), ylim=c(mn,mx),
pch=16, cex=0.5, col="gray10")
curve(1*x, from=mn, to=mx, add=T, lty=2, lwd=2, col="blue")
bp1 <- round(mean(rep > obs),2)
loc <- ifelse(bp1 < 0.5, "topleft", "bottomright")
legend(loc, legend=bquote(p[B]==.(bp1)), bty="n", cex=2)
if (ratio==TRUE){
t.rep <- samps["tvm.rep"][[1]][ind,]
t.obs <- samps["tvm.obs"][[1]][ind,]
# plot variance/mean ratio
hist(t.rep, nclass=50,
xlab="variance/mean ", main=NA, axes=FALSE)
abline(v=t.obs, col="red")
axis(1); axis(2)
}
return(list('Bayesian p-value'=bp1))
}
# Counts of adults (nonbreeders+breeders)
fit.check(out, ratio=F,
name.rep="dmape.rep",
name.obs="dmape.obs",
ind=1,
lab="Adults(Breeder+Nonbreeder)- Poisson", jit=300)
## $`Bayesian p-value`
## [1] 0.66# Counts of first-year fledglings, ind=2
# Poisson failed fit test bp=0
# So we assigned negative binomial only for Punta Cana.
# Los Haitises excluded from neg bin because no zeroes in counts.
fit.check(out, ratio=F,
name.rep="dmape.rep",
name.obs="dmape.obs",
ind=2,
lab="First-year counts\nNeg binomial-Poisson", jit=300)
## $`Bayesian p-value`
## [1] 0.57# Productivity
fit.check(out, ratio=F,
name.rep="f.dmape.rep",
name.obs="f.dmape.obs",
ind=1,
lab="Productivity-Neg binomial", jit=300)
## $`Bayesian p-value`
## [1] 0.833.2 Examine Traceplots, Compare Posteriors with Priors
Traceplots provide evidence that models adequately converged.
MCMCtrace(post2, pdf=FALSE, params= "mus", Rhat=TRUE, priors=runif(20000, 0, 1), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "betas", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "deltas", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "gamma", Rhat=TRUE, priors=runif(20000, -20, 20), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "sds2", Rhat=TRUE, priors=rexp(20000, 1), post_zm=FALSE)
MCMCtrace(post2, pdf=FALSE, params= "R2")
